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FIRE-PLUME-GENERATED  CEILING  JET  CHARACTERISTICS  AND  CONVECTIVE 
HEAT  TRANSFER  TO  CEILING  AND  WALL  SURFACES  IN  A TWO-LAYER  ZONE- 
TYPE  FIRE  ENVIRONMENT:  UNIFORM  TEMPERATURE  CEILING  AND  WALLS 

Leonard  Y.  Cooper 


ABSTRACT 

It  has  been  determined  by  Sandia  National  Laboratories  and  the  US  Nuclear  Regulatory  Commission 
that  the  use  of  deterministic,  multi-room,  zone-type  fire  modeling  technology  could  enhance  the 
reliability  of  their  recent  reactor  safety  risk  studies.  These  studies  are  confined  to  the  relatively  early 
detection  times  of  fire  development  when  fire-driven  ceiling  jets  and  gas-to-ceiling  convective  heat 
transfer  are  expected  to  play  a particularly  important  role  in  room-to-room  smoke  spread  and  in  the 
response  of  near-ceiling  mounted  detection  hardware.  A parameter  of  concern  in  these  risk  analyses 
is  the  location  of  the  fire  within  the  space  of  fire  origin.  One  goal  of  the  analyses  is  to  determine 
the  significance  to  risk  of  this  fire-position  parameter. 

This  work  presents  a model  to  predict  the  instantaneous  rate  of  convective  heat  transfer  from  fire 
plume  gases  to  the  overhead  ceiling  surface  in  a room  of  fire  origin.  The  room  is  assumed  to  be  a 
rectangular  parallelopiped  and,  at  times  of  interest,  ceiling  temperatures  are  simulated  as  being 
uniform.  Also  presented  is  an  estimate  of  the  convective  heat  transfer,  due  to  ceiling-jet-driven  wall 
flows,  to  both  the  upper  and  lower  portions  of  the  walls.  The  effect  on  the  heat  transfer  of  the 
location  of  the  fire  within  the  room  is  taken  into  account.  Finally  presented  is  a model  of  the  velocity 
and  temperature  distributions  in  the  ceiling  jet. 

The  model  equations  were  used  to  develop  an  algorithm  and  associated  modular  computer  subroutine 
to  carry  out  the  indicated  heat  transfer  calculations.  The  subroutine  is  written  in  FORTRAN  77  and 
called  CEILHT.  The  algorithm  and  subroutine  are  suitable  for  use  in  two-layer  zone-type 
compartment  fire  model  computer  codes.  The  subroutine  was  tested  for  a variety  of  fire 
environments  involving  a 107W  fire  in  a 8m  x 8m  x 4m  high  enclosure.  While  the  calculated  results 
were  plausible,  it  is  important  to  point  out  that  CEILHT  simulations  have  not  been  experimentally 
validated. 


Keywords:  building  fires;  ceiling  jets;  compartment  fires;  computer  models;  fire  models;  heat 

transfer;  mathematical  models;  zone  models 


Introduction 


It  has  been  determined  by  Sandia  National  Laboratories  and  the  US  Nuclear  Regulatory  Commission 
that  the  use  of  deterministic,  multi-room,  zone-type  fire  modeling  technology  could  enhance  the 
reliability  of  their  recent  reactor  safety  risk  studies.  These  studies  are  confined  to  the  relatively  early 
detection  times  of  fire  development  when  fire-driven  ceiling  jets  and  gas-to-ceiling  convective  heat 
transfer  are  expected  to  play  a particularly  important  role  in  room-to-room  smoke  spread  and  in  the 
response  of  near-ceiling  mounted  detection  hardware.  A parameter  of  concern  in  these  risk  analyses 
is  the  location  of  the  fire  within  the  space  of  fire  origin.  One  goal  of  the  analyses  is  to  determine 
the  significance  to  risk  of  this  fire-position  parameter.  The  purpose  of  this  work  is  to  develop  and 
present  a means  of  simulating  the  heat  transfer  and  flow  phenomena  so  that  they  can  be  included 
in  two-layer  zone-type  fire  model  computer  codes. 

Consider  a room  of  fire  origin  in  a multi-room  facility.  At  some  arbitrary  instant  during  a dynamic 
fire  scenario  assume  that  the  fire-generated  environment  in  the  room  can  be  reasonably  simulated 
by  upper  and  lower  uniform-property  layers  where  the  two  layers  are  stably-stratified  (i.e.,  the  upper 
layer  temperature  is  greater  than  the  lower  layer  temperature). 

Consistent  with  several  two-layer  zone-type  compartment  fire  model  computer  codes  (see  e.g., 
References  [1],  [2],  and  [3])  assume  that  it  is  reasonable  to  simulate  the  room’s  ceiling  temperature 
as  uniform.  One  objective  of  this  work  is  to  develop  an  algorithm  and  associated  computer 
subroutine,  useable  in  such  fire  models,  to  predict  the  instantaneous  rate  of  convective  heat  transfer 
from  the  fire  plume  to  the  ceiling  surface.  In  developing  the  heat  transfer  simulation  it  is  desired  to 
account  for  the  location  of  the  fire  within  the  room.  Estimates  of  the  rates  of  convective  heat 
transfer  to  the  upper-layer  portion  of  the  wall  surfaces  (i.e.,  the  portion  submerged  by  the  upper 
layer)  and  lower-layer  portion  of  the  wall  surfaces  will  also  be  provided.  A second  objective  of  this 
work  is  to  provide  an  estimate  of  the  velocity  and  temperature  distributions  in  the  ceiling  jet  which 
is  formed  by  plume/ceiling  impingement  and  radial  spread  of  the  plume  gases. 

Let  the  room  be  a rectangular  parallelopiped  where  the  geometry  is  specified  and,  at  a time  of 
interest,  assume  that  the  fire  environment  in  the  room  and  the  characteristics  of  the  fire  itself  have 
been  determined. 


Room  Geometry,  Characterization  of  the  Instantaneous  State  of  the  Fire,  the  Fire  Plume,  and  the 
Fire-Generated  Environment 

Define  a right-handed  cartesian  coordinate  system  with  axes  X,Y,Z  and  with  origin  at  a lower  corner 
of  the  room,  i.e.,  at  the  floor  elevation.  Orient  the  axes  so  that  the  Z axis  is  directed  upward  and 
the  X and  Y axes  are  aligned  along  the  junction  of  the  walls  and  the  floor.  Designate  the  co- 
ordinates of  corner  of  the  room  furthest  from  the  origin  of  the  axes  as  XWALL,  YWALL,  and  ZCEIL. 
These  are  the  length,  width,  and  ceiling  height  of  the  room,  respectively.  Designate  the  co-ordinates 
of  the  center  of  the  base  of  the  fire  as  XnRE,  YFIRE,  and  ZFIRE. 

Characterize  and  specify  the  state  of  the  fire-generated  environment  in  the  room  by:  the  average 
densities  and  assumed  uniform  temperatures  of  the  upper  and  lower  gas  layers,  py,  Ty,  andpL,  TL, 
respectively;  the  temperature  of  the  ceiling  surface  and  the  average  temperature  of  the  wall  surfaces. 
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Tceil  and  Twat  t , respectively;  and  the  elevation  above  the  floor  of  the  interface  between  the  two 
gas  layers,  ZEAYER.  A P^an  yiew  and  an  elevation  view  of  the  geometry  is  depicted  in  Figure  1. 

Define  and  specify  QpcoNV  as  t^ie  portion  of  the  total  fire  energy  release  rate  convected  in  the  fire 
plume  (i.e.,  the  total  rate  of  combustion  energy  release  less  the  rate  of  energy  radiated  from  the 
combustion  zone  and  the  plume).  If  the  fire  is  below  the  layer  interface,  ZnRE  < ZEAYER,  then 
characterize  the  mass  flow  rate  and  the  convected  enthalpy  flow  rate  in  the  plume  at  ZLAYER  as 
Mplume  and  QpLUME’  respectively. 

Let  Mplume  be  specified.  Then,  assuming  that  the  fuel  flow  rate  is  negligible  compared  to  Mplume, 
Qplume  = MPLUMECpTL  + Qfconv  (!) 


Overview  of  Plume/Upper-Layer  Interaction  and  Ceiling  Heat  Transfer 

The  buoyant  fire  plume  rises  from  ZEIRE  toward  the  ceiling.  When  the  fire  is  below  the  layer 
interface,  its  mass  and  enthalpy  flow,  Mplume  and  QpLUME’  are  assumed  to  be  deposited  into  the 
upper  layer  at  ZLAYER.  Having  penetrated  the  interface,  a portion  of  the  plume  typically  continues 
to  rise  toward  the  ceiling.  As  it  impinges  on  the  ceiling  surface,  the  plume  gases  turn  and  form  a 
relatively  high  temperature,  high  velocity,  turbulent  ceiling  jet  which  flows  radially  outward  along  the 
ceiling  and  transfers  heat  to  the  relatively  cool  ceiling  surface.  The  ceiling  jet  is  cooled  by  convection 
and  the  ceiling  material  is  heated  indepth  by  conduction.  The  convective  heat  transfer  rate  is  a 
strong  function  of  the  radial  distance  from  the  point  of  plume/ceiling  impingement,  reducing  rapidly 
with  increasing  radius.  It  is  dependent  also  on  the  characteristics  of  the  plume  immediately  upstream 
of  ceiling  impingement. 

The  relatively  high  temperature  ceiling  jet  is  blocked  eventually  by  the  relatively  cool  wall  surfaces. 
IN  THIS  WORK  IT  IS  ASSUMED  THAT  THE  DISTANCE  BETWEEN  THE  PLUME/CEILING 
IMPINGEMENT  POINT  AND  THE  CLOSEST  WALL  SURFACE  IS  NEVER  LESS  THAN  0.2(ZCEIL 
- ZFIRE).  Ceiling-jet/wall-surface  impingement  leads  to  stagnation-point-type  enhancement  of  jet-to- 
surface  heat  transfer  [4].  The  ceiling  jet  then  turns  downward  and  outward  in  a complicated  flow 
along  the  vertical  wall  surfaces  (see  [5]  and  [6]).  Convective  wall-flow-to-wall-surface  heat  transfer 
continues.  The  descent  of  the  wall  flows  and  the  heat  transfer  from  them  are  stopped  eventually  by 
upward  buoyant  forces.  They  are  then  buoyed  back  upward  and  mix  finally  with  the  upper  layer. 

Let  qEE1E(X,Y)  be  the  instantaneous  flux  of  convective  heat  transfer  from  the  ceiling  jet  gases  to  the 
ceiling  surface  and  let  flcnL  ave  and  Qceil  be  its  average  value  and  the  value  of  its  integral, 
respectively,  over  the  entire  ceiling  surface.  The  first  objective  of  this  work  is  to  determine 


QcEIL  “ (XWALLYWALL)(icEIL,AVE  “ 


XWALL  YWALL 

/ / qcEIL(X,Y)dXdY 

0 0 


(2) 
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Properties  of  the  Plume  in  the  Upper  Layer  When  ZFIRE  < Zj^yer 

Consider  times  when  the  elevation  of  the  fire  is  below  the  interface,  i.e.,  when  ZFIRE  < ZLAYER. 

As  the  plume  flow  enters  the  upper  layer  the  forces  of  buoyancy  which  drove  the  plume  toward  the 
layer  interface  (i.e.,  as  a result  of  relatively  high-temperature,  low-density  plume  gases  being 
submerged  in  a relatively  cool,  high-density  lower-layer  environment)  are  reduced  immediately 
because  of  the  increase  of  temperature  of  the  upper  layer  environment  over  that  of  the  lower  layer. 
The  continued  ascent  of  the  plume  gases  will  be  less  vigorous,  i.e.,  at  reduced  velocity,  and  of  higher 
temperature  than  it  would  have  been  in  the  absence  of  the  layer.  Indeed,  some  of  the  penetrating 
plume  flow,  will  actually  be  at  a lower  temperature  than  Ty.  The  upper  layer  buoyant  forces  on  this 
latter  portion  of  the  flow  will  actually  retard  and  possibly  stop  its  subsequent  rise  to  the  ceiling. 

The  simple  point-source  Gaussian-distribution  plume  model  of  [7]  will  be  used  to  simulate  the  plume 
flow,  first  immediately  below,  or  upstream  of  the  interface,  and  then  throughout  the  depth  of  the 
upper  layer  itself. 

Consider  a Reference-[7]-type  point  source  of  buoyancy,  with  elevation  below  ZLAYER.  Then  the 
plume  above  such  a source  will  be  taken  to  be  equivalent  to  the  plume  of  our  fire  (in  the  sense  of 
having  identical  mass  and  enthalpy  flow  rates  at  the  interface)  if  the  point  source  strength  is 
Qfconv  ar*d  the  elevation  of  the  equivalent  source,  ZEQ,  satisfies 

^PLUME  = 0-21PlS1;/2(ZlayER  ' ZEq)5/2QeQ1/3  (3) 

where  QEq,  a dimensionless  measure  of  the  strength  of  the  fire  plume  at  ZLAYER,  is  defined  by 

QeQ  = QFCONv/l/)LCpTLg1/"(ZLAYER  ‘ ZEq)5/“]  (4) 

and  where  Cp  is  the  specific  heat  at  constant  pressure  of  the  entrained  gas,  which  can  be  assumed 
to  be  identical  to  that  of  air. 

Solving  (3)  and  (4)  for  ZEq  and  Qeq  in  terms  of  known  variables  leads  to 


Qeq  - [^^Qfconv^^p^l^plume)]3^  (5) 

ZEQ  = Z LAYER  ' [QfCONv/CQeO^L^p^lS1^)]^5  (6) 

As  the  plume  crosses  the  interface,  the  fraction,  M*,  of  Mplume  which  is  still  buoyant  relative  to  the 
upper  layer  environment  and  presumably  continues  to  rise  to  the  ceiling,  entraining  upper  layer  gases 
along  the  way,  is  predicted  in  [8]  to  be 
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f° 

M*  = 

l (1.04599o-  + 0.36039  1c72)/(1.  + 1.37748a  + 0.360391a2) 


if  -1  < a < 0; 

(7) 

if  a > 0 


where  the  dimensionless  parameter  a is  defined  as 

a — (1  - a + CxQ'qm)/(«  - 1)  (8) 

a = TuATl  ; Cx  = 9.115  (9) 


If  M*  > 0,  Reference  [8]  further  identifies  the  parameters  necessary  to  describe  plume  flow  continua- 
tion in  the  upper  layer  (i.e.,  between  ZLAYER  and  ZCEIL)  according  to  a Reference-[7]-type  point 
source  plume.  It  has  been  determined  there  that  this  plume  can  be  modeled  as  being  driven  by  a 
non-radiating  continuation  buoyant  point  source  of  strength  Q’  located  a distance 


H “ ZCEIL  ' ^SOURCE  > ZCEIL  * ZFIRE  (10) 

below  the  ceiling  in  a (downward-)  extended  upper  layer  environment  of  temperature  Ty  and  density 
Pjj.  The  relevant  parameters  predicted  in  [8]  are 


Q - QpyoNv^M '7(1  + a)  (11) 

^SOURCE  = -^LAYER  - (ZLAYer  - ZEQ)a3/5M  ' "/5[(1  + a)/a]1/:i  (12) 

The  fire  and  the  equivalent  source  in  the  lower  layer  and  the  continuation  source  in  the  upper  layer 
are  depicted  in  Figure  2.  Times  during  a fire  simulation  when  Eq.  (8)  predicts  a > > 1 are  related 
to  states  of  the  fire  environment  when  the  temperatures  above  TL  of  the  radial  Gaussian  temperature 
distribution  of  the  plume  flow,  at  the  elevation  of  interface  penetration,  is  predicted  to  be  mostly 
much  larger  than  (Ty  - TL).  Under  such  circumstances  the  penetrating  plume  flow  is  still  very 
strongly  buoyant  as  it  enters  the  upper  layer.  The  plume  will  continue  to  rise  to  the  ceiling  and  to 
drive  ceiling  jet  convective  heat  transfer  at  rates  which  differ  only  slightly  (on  account  of  the  elevated 
temperature  upper  layer  environment)  from  the  heat  transfer  rates  which  could  occur  in  the  absence 
of  an  upper  layer. 

Conditions  where  Eq.  (8)  predicts  a < 0 are  related  to  times  during  a fire  scenario  when  the 
temperature  of  the  plume  at  the  elevation  of  interface  penetration  is  predicted  to  be  uniformly  less 
than  Ty.  Under  such  circumstances  the  penetrating  plume  flow  is  nowhere  positively  (i.e.,  upward) 
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buoyant  as  it  enters  the  upper  layer.  Thus,  while  all  of  this  flow  is  assumed  to  enter  and  mix  with 
the  upper  layer,  it  is  assumed  further  that  none  of  it  rises  to  the  ceiling  in  a coherent  plume,  i.e.,  Q’ 
= 0.  For  this  reason,  when  a < 0 the  model  precludes  to  existence  of  any  significant  ceiling  jet  flow 
and  associated  convective  heat  transfer  to  the  ceiling  surface. 

The  above  analysis  assumes  that  ZnRE  < ZLAYER.  However,  at  the  onset  of  a fire  scenario  it  is 
typical  that  ZFIRE  < ZLAYER  = ZCEIL  and  a , a and  M*  of  Eqs.  (7)-(9),  which  depend  on  the 
indeterminate  initial  value  of  Tu?  are  themselves  undefined.  The  situation  when  Zj  AYER  = ZCEIL 
is  properly  taken  into  account  if  Q’  = QEconv  anc*  zsource  = zeq- 

General  Properties  of  the  Plume  in  the  Upper  Layer  and/or  Near  the  Ceiling 

When  ZnRE  < ZLAYER  < ZCEIL,  the  results  of  Eqs.  (11)  and  (12)  allow  a description  of  the  fire- 
driven  plume  dynamics  in  the  upper  layer,  up  to  plume/ceiling  impingement  at  ZLAYER,  i.e.,  by  using 
the  continuation  point  source  and  the  plume  model  of  [7]  in  the  extended  upper  layer.  When 
z layer  < zfire>  the  same  plume  model  is  used,  but  Q’  and  Z^qurce  are  taken  to  be  the  convected 
buoyant  strength  of  the  fire  and  the  actual  fire  elevation,  respectively.  Finally,  when  ZLAYER  = 
ZCEil>  the  plume  dynamics  near  the  ceiling  is  described  by  the  equivalent  point  source,  which  is  in 
the  lower  layer,  and  the  point  source  plume  model. 

All  cases  can  be  treated  with  the  following  final  versions  of  Eqs.  (11)  and  (12) 


Q’  = 


| QfcONV0^ ' /(l  + a) 
l Qfconv 


ZFIRE  < Z LAYER  < ZCEIL 
^ ZFIRE  - Z LAYER  or  Z LAYER  = ZCEIL 


(IF) 


ZSOURCE  “ 


ZEQ  ‘f  ZLAYER  “ ZCEILi 

Z LAYER  " (ZLAYER  ' ZEQ)o:3/5M  I'-/5[(l  + a)/c r]1/5 

. 

if  ZFIRE  < Z LAYER  < ZCEIL’ 
. ZFIRE  if  ZLAYER  - ZFIRE  < ZCEIL 


(12’) 


where  ZEq,  M*,  a,  and  a are  calculated  from  Eqs.  (5)-(9). 

Calculating  the  Convective  Heat  Transfer  to  the  Ceiling 

When  the  fire  is  below  the  interface  and  the  interface  is  below  the  ceiling  the  method  of  [9]  will  be 
used  to  calculate  qEER.  The  method  of  [9]  was  developed  to  treat  generic,  confined-ceiling,  room 
fire  scenarios.  The  confined  ceiling  problem  is  solved  by  applying  the  unconfined  ceiling  heat  transfer 
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solution  to  the  problem  of  an  Eq.  (1 1’)-(12’)  equivalent  upper  layer  source  in  an  extended  upper  layer 
environment  [10,  11].  When  the  fire  is  above  the  interface,  the  unconfined  ceiling  methodology 
applies  directly.  This  solution  technique  has  been  used  previously  in  [12]  and  [13]. 

The  equivalent,  extended-upper-layer,  unconfined-ceiling  flow  and  heat  transfer  problem  is  depicted 
in  the  right-hand  sketch  of  Figure  2.  It  involves  the  equivalent  Q’  heat  source  of  Eq.  (IF)  located 
a distance  H below  the  ceiling  surface  in  an  extended  lower  layer  environment  where  H is  found  from 
Eqs.  (10)  and  (12’).  The  objective  here  is  to  estimate  the  q£EIL(X,Y)  distribution  for  the  scenario 
in  the  latter  sketch,  and  then  its  Eq.  (2)  integral  over  the  ceiling  surface. 

From  [10]  and  [11] 


9cEIl(X’Y)  ~ hL(TAD  " TCEIl) 


(13) 


where  Tad,  a characteristic  ceiling  jet  temperature,  is  the  temperature  that  would  be  measured 
adjacent  to  an  adiabatic  lower  ceiling  surface,  and  hE  is  a heat  transfer  coefficient.  hE  and  T/VD  are 
given  by 


j 8.82ReH-1/2Pr'2/3[l  - (5.0  - 0.284ReHa2)(r/H)]  if  0 < r/H  < 0.2; 

hL/h  = (14) 

l 0.283ReH'0  3Pr'2/3(r/H)'1'2(r/H  - 0.0771)/(r/H  + 0.279)  if  0.2  < r/H 


(T^  - TuMTuQif3)  = 


[ 10.22  - 14.9r/H 
{ 8.39f(r/H) 


if  0 < r/H  < 0.2; 


(15) 

if  0.2  < r/H 


where 


f(r/H)  = [1  - 1.10(r/H)°'8  + 0.808(r/H)L6]/ 

[1  - 1.10(r/H)°'8  + 2.20(r/H)1'6  + 0.690(r/H)14] 
r = [(X  - XnRE)2  + (Y  - YF1RE)2]1/2 

fi  = PuCpS1/2h1/2°h1/3;  ReH  = g1/2H3/2QA1/3/vu;  Q£  = Q7buCpTu(gH)1/2H2] 


(16) 

(17) 

(18) 


In  the  above,  Pr  is  the  Prandtl  number  (taken  to  be  0.7)  and  is  the  kinematic  viscosity  of  the 
upper  layer  gas  which  is  assumed  to  have  the  properties  of  air.  Also,  Qj1),  a dimensionless  number, 
is  a measure  of  the  strength  of  the  plume  and  ReH  is  a characteristic  Reynolds  number  of  the  plume 
at  the  elevation  of  the  ceiling. 
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The  following  estimate  for  Vy  of  air  [14]  will  be  used  when  computing  ReH  from  Eqs.  (18) 


vv  = [0.04128(10‘7)Tu5/2]/[Tu  + 110.4]  vu  in  m2/s,  Ty  in  K (19) 

Eqs.  ( 13)-(  19)  are  valid  when  ZLAYER  < ZCEIL  and  Ty  is  well  defined.  When  ZEAYER  = ZCEIL  and 
Ty  is  undefined,  the  results  are  also  valid  if  Ty  is  set  equal  to  TL  in  these  equations.  This  will  yield 
the  correct  limiting  result  for  the  convective  heat  transfer  to  the  ceiling,  namely,  convective  heat 
transfer  to  the  ceiling  from  an  unconfined  ceiling  jet  in  a lower  layer  environment. 

Integrating  qEEIL(X,Y)  according  to  Eq.  (2)  yields  the  desired  results  for  QCEIL  and  qEEIL  AVE. 


Calculating  the  Average  Convective  Heat  Transfer  to  the  Upper  and  Lower  Wall  Surfaces 

As  mentioned  above  Eq.  (2),  the  ceiling  jet  is  blocked  eventually  by  the  relatively  cool  wall  surfaces 
where  ceiling-jet/wall-surface  impingement  leads  to  a stagnation-point  type  of  heat  transfer 
enhancement  there.  In  [4]  it  is  estimated  that  for  fires  relatively  close  to  a wall  (i.e.,  a wall  near  to, 
but  outside,  the  presently  restricted  r/H  = 0.2  circle,  where  ceiling  jet  velocities  are  at  or  near  to  their 
largest  amplitudes)  the  rate  of  heat  transfer  to  the  wall,  q^ALL  is  more  than  two  times  that  of  the 
jet-to-local-ceiling-surface  heat  transfer. 

Locations  of  Normal  Ceiling-, Tet/Wall  Impingement.  For  the  fire  scenarios  being  studied  here,  there 
are  four  locations  where  ceiling-jet/wall  impingement  involves  ceiling  jet  velocities  normal  to  the  wall. 
These  locations  are  identified  in  Figure  1 as  STP1,  STP2,  STP3,  and  STP4.  At  such  stagnation 
points,  i.e.,  at  the  ceiling/wall  junction,  an  estimate  of  the  heat  transfer  to  the  wall  is  [4] 


9wall,st  ~ ^st(^ad  " Twall)  (-0) 

where 

hST/h  = 0.94ReH-a42/[(rST/H)Pr]  (21) 

rST  is  the  value  of  r of  Eq.  (16)  at  the  stagnation  point  (e.g.,  rST  = YEIRE  for  STP1  in  Figure  1),  and 
Taj},  R,  and  ReH  are  given  in  Eqs.  (15)  and  (17). 

In  [5]  it  is  estimated  that  at  locations  of  normal  jet/wall  impingement  the  penetration  depth  of  the 
downward-directed  wall  flow  is  approximately  0.8H  or  ZCEjE  (i.e.,  penetration  to  the  floor),  whichever 
is  smaller,  and  that  this  result  is  relatively  invariant  with  the  value  of  rST/H.  Here  it  is  assumed  that 
the  heat  transfer  goes  to  zero  at  this  penetration  depth,  and  remains  zero  below  this  elevation.  It 
is  assumed  further,  that  the  heat  transfer  to  the  wall  is  linear  with  elevation,  going  from  4\vall,ST 
at  the  ceiling  to  zero  at  Z = ZCEIE  - 0.8H.  (It  is  possible  for  the  latter  elevation  to  be  negative,  i.e., 
below  the  floor,  because  of  the  fact  that  when  the  fire  is  in  the  lower  layer,  H is  based  on  the 
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elevation  of  an  equivalent  point  source  of  buoyancy  in  the  extended  upper  layer.  Consistant  with 
the  inequality  of  Eq.  (10),  this  elevation  can  actually  be  below  the  floor  elevation,  Z = 0.  When  such 
is  the  case,  the  assumed  linear  distribution  for  wall  heat  transfer  is  only  used  between  the  floor  and 
ceiling  elevations. 

The  above  paragraph  leads  to  the  following  estimate  for  the  average  convective  heat  transfer  on  the 
upper  wall  segment  (between  ZCEIL  and  ZLAYER),  Owaix  NORM  U > an^  on  the  lower  wall  segment 
(between  ZLAYER  and  the  floor),  Awall  norm  L’  along  vertical  lines  on  a wall  where  normal  ceiling- 
jet/wall  impingement  occurs: 


0.8H  < ZCEIL  - ZLAYER: 

fl WALL, NORM, U = (fl\VALL,ST/2)(ZCEIL  * O.SH)/(ZCEIL  ' ZLAYEr) 
fl  WALL, NORM, L = 0 

ZCEIL  ' Z LAYER  - < ZCEIL: 


(22) 

(23) 


Z 


WALL, NORM, U 

WALL, NORM, L 

ceil  — 0-8H: 

fl  WALL, NORM, U 


AwalestC1  "(zceil  " zlayer)/(1-6H)] 


(flWALL.SI72)  [0-8H  - (ZCEIL  - ZLAYER)]/ZLAYER 


flWALL.STt1  '(ZCEIL  " ZLAYErV(1-6H)] 


(24) 

(25) 


(26) 


fl WALL, NORM, L “ (AwaLL,St/2)[1-6H  * (ZCEIL  * ZLAYEr)  ’ ZCEIlJ/(®-^)  (27) 

Strategy  for  Estimating  Heat  Transfer  to  the  Walls.  As  noted,  Eqs.  (22)-(27)  can  be  applied  to 
vertical  lines  passing  through  points  STP1,  STP2,  STP3,  and  STP4.  The  following  strategy  is  adopted 
for  estimating  the  average  rate  of  ceiling-jet-driven  convective  heat  transfer  to  the  upper  and  lower 
layer  wall  segments  of  the  room: 

1.  Use  Eqs.  (20)  and  (21)  to  calculate  q\YAEE  at  the  top  of  the  walls  near  the  points  STP1,  STP2, 
STP3,  and  STP4.  Then  use  Eqs.  (22)-(27)  to  calculate  the  average  q\vALL  al°nS  the  upper-  and 
lower-layer  segments  of  vertical  lines  through  these  points. 

2.  Estimate  q(vAEL  at  the  top  of  the  walls  near  the  four  room  corners  of  the  room  by  assuming 
normal  ceiling-jet/wall  impingement  there  and  using  Eqs.  (20)  and  (21).  For  example,  use  these 


9 


equations  with  rST  = (XpIRE  4-  YpIRE)1/2  to  estimate  the  heat  ransfer  near  the  top  of  the  walls 
at  X = Y = 0.  Then  use  Eqs.  (22)-(27)  to  calculate  the  average  q^ALL  al°ng  the  upper-  and 
lower-layer  segments  of  the  vertical  lines  which  define  these  corners. 

3.  For  an  upper  or  lower  wall  segment  between  STP1,  STP2,  STP3,  or  STP4,  and  an  adjacent 
corner,  estimate  the  average  value  of  q^ALL  as  the  average  of  the  previously  computed  average 
values  along  the  vertical  bounding  "end-lines."  Calculate  this  average  for  each  of  the  eight 
upper  segments  and  eight  lower  segments. 

4.  For  each  of  the  upper-layer  segments  use  the  area  of  the  segment  and  the  previously  calculated 
average  q\vALE  to  calculate  the  eight  contributions  to  the  total  rate  of  upper-layer  wall  heat 
transfer.  Sum  these  contributions  and  obtain  finally  4walluavE’  the  average  tlux  of  heat 
transfer  to  the  upper-layer  portions  of  the  walls.  Carry  out  analogous  calculations  for 
4wall  l AVE’  the  average  flux  of  heat  transfer  to  the  lower-layer  portions  of  the  walls. 

5.  Add  the  rates  of  heat  transfer  to  all  sixteen  upper  and  lower  wall  segments  and  obtain,  QwalL’ 
the  total  rate  of  heat  transfer  to  the  wall  surfaces. 

The  details  of  the  above  algorithm  are  presented  in  Appendix  A. 

It  is  important  to  point  out  that  in  compartment  fires  there  are  fluid  dynamic  effects  other  than 

ceiling-jet/wall  impingement  that  can  lead  to,  and  dominate  the  convective  heat  transfer  to  walls  (see, 

e.g.  [6]  and  [15]).  These  effects  are  not  considered  in  this  work. 


The  Computer  Subroutine  CEILIIT  for  Carrying  out  the  Heat  Transfer  Calculations 

The  model  equations  of  the  last  two  sections  were  used  to  develop  an  algorithm  and  associated 
modular  computer  subroutine  to  carry  out  the  indicated  heat  transfer  calculations.  The  subroutine 
is  written  in  FORTRAN  77  and  is  called  CEILHT.  A summary  flow  diagram  for  the  algorithm  is 
presented  in  Figure  3.  The  algorithm  and  subroutine  are  suitable  for  use  in  two-layer  zone-type 
compartment  fire  model  computer  codes. 

CEILHT,  which  uses  the  two-dimensional  numerical  integration  subroutine  ADAPT  [16],  is  presented 
in  Appendix  B.  CEILHT  has  been  tested  for  a variety  of  instantaneous  fire  environments  involving 
a 107W  fire  in  a 8m  x 8m  x 4m  high  enclosure.  While  the  calculated  results  were  plausible,  it  is 
important  to  point  out  that  CEILHT  simulations  have  not  been  experimentally  validated. 


Velocity  and  Temperature  Distributions  in  the  Ceiling  Jet 

It  is  an  objective  of  this  work  to  simulate  the  velocity  and  temperature  distributions  in  the  ceiling  jet. 
The  major  interest  is  to  develop  a general  capability  for  predicting  the  dynamic  environment  local  to 
near-ceiling-deployed  fire  safety  devices,  e.g.  fusible  sprinkler  links.  Knowledge  of  this  environment 
would  be  used  in  predictions  of  their  thermal  response  to  the  fire-generated  environment  (see,  e.g., 
[13]  and  [17]). 
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For  relatively  smooth-ceiling  configurations,  assumed  to  be  representative  of  the  rooms  of  fire  origin 
studied  in  this  work,  the  ceiling  jet  flows  outward  radially  from  the  point  of  plume/ceiling 
impingement  and  its  gas  velocity  and  temperature  distributions,  Va  and  TCJ,  respectively,  are  a 
function  of  radius  from  the  impingement  point,  r,  distance  below  the  ceiling,  z,  and  time.  Note  that 
it  has  been  shown  with  some  generality  that  plume-driven  ceiling  jets  are  inertially  dominated  flows 
from  small  to  moderate  r/H  values  [18]. 

The  Velocity  Distribution.  Let  z be  the  distance  downward  from  the  ceiling  surface,  i.e.,  if  Z is  the 
elevation  above  the  floor, 


z “ ^CEIL  - Z (28) 

Define  the  plume/ceiling  impingement  stagnation  zone  by  r/H  < 0.2.  Outside  the  stagnation  zone 
and  at  a given  r,  Vqj  rises  rapidly  from  zero  at  the  ceiling’s  lower  surface,  z = 0,  to  a maximum, 
VMAX,  at  a distance  z = 0.235,  5(r)  being  the  distance  below  the  ceiling  where  V/VMAX  = 1/2  [5]. 
In  this  region  outside  the  stagnation  zone,  VCJ  can  be  estimated  from  [5] 


When  r/H  > 0.2: 


VCj/VMAX  = 


(8/7)[z/(0.235)]1/7{  1 - [z/(0.235)]/8}  if  0 < z/(0.235)  < 1; 

(29) 

cosh'2{(0.23/0.77)arccosh(21/2)[z/(0.235)  - 1]}  if  1 < z/(0.235) 


VMAX/V  = 0.85(r/H)-11;  5/H  = 0.10(r/H)a9;  V = g1/2H1/2Q£1/3  (30) 

where  is  defined  in  Eq.  (18).  Vqj/Vj^^  per  Eq.  (29)  is  plotted  in  Figure  4. 

Inside  the  stagnation  zone,  the  fire-driven  flow  is  changing  directions  from  an  upward-directed  plume 
flow  to  a outward-directed  ceiling-jet-type  flow.  There  the  flow  velocity  involves  generally  a 
significant  vertical  as  well  as  radial  component  of  velocity  with  the  characteristic  velocity  being  the 
maximum  ceiling  jet  velocity  at  r/H  = 0.2 


When  0 < r/H  < 0.2: 


I vcj  I = °rder  of  = °-2)] 


(31) 


Other  considerations  are  required  for  a more  accurate  description  of  the  velocity  in  the  stangnation 
zone. 
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The  Temperature  Distribution.  Outside  of  the  stagnation  zone,  i.e.,  where  r/H  > 0.2,  and  at  a given 
value  of  r,  TCJ  changes  very  rapidly  from  the  temperature  of  the  ceiling’s  lower  surface  TCEI[ , at  z 
= 0,  to  a local  maximum,  TMAX,  somewhat  below  the  ceiling  surface.  We  assume  that  this  maximum 
value  of  TCJ  occurs  at  the  identical  distance  below  the  ceiling  as  does  the  maximum  of  VCJ,  i.e.,  at 
z = 0.235.  Below  this  elevation,  TCJ  drops  with  increasing  distance  from  the  ceiling  until  it  reaches 
the  upper  layer  temperature,  Ty.  In  this  latter  outer  region  of  the  ceiling  jet,  the  shape  of  the 
normalized  TCJ  distribution,  (Tqj  - Ty)/(TMAX  - Ty),  will  have  the  same  characteristics  as  does  that 
of  Va/VMAX.  Also,  since  we  are  dealing  with  a turbulent  boundary  flow  it  is  reasonable  to  expect 
that  the  characteristic  thicknesses  of  the  outer  region  of  both  the  velocity  and  temperature 
distributions  are  the  same,  being  dictated  there  by  the  distribution  of  the  turbulent  eddies. 

For  the  above  reasons  we  approximate  the  velocity  and  temperature  distributions  as  being  identical 
in  the  outer  region  of  the  ceiling  jet  flow,  0.235  < z.  In  the  inner  region  of  the  flow,  between  z = 
0 and  0.235,  we  approximate  the  normalized  temperature  distribution  by  a quadratic  function  of 
z/(0.235),  requiring  TCJ  = TCE1L  at  z = 0 and  Tqj  = TMAX,  dTCJ/5z  = 0 at  z = 0.235.  Thus, 


When  r/H  > 0.2: 


0S  + 2(1  - 0s)[z/(O.235)]  - (1  - 0s)[z/(O.235)]2 


® - (TCj-Tu)/ (Tmax'^u)  - 


if  0 < z/(0.235)  < 1 (32) 


. vo//vmax 


if  1 < z/(0.235) 


0S  - ©(Tyj  - Tceil)  - (TCEIL  - Ty)/(TMAX  - Ty)  (33) 

Note  that  0S  will  be  negative  where  the  ceiling  surface  temperature  is  less  than  the  upper  layer 
temperature,  for  example,  relatively  early  in  the  fire  when  the  original  ambient-temperature  ceiling 
surface  has  not  yet  reached  the  average  temperature  of  the  growing  upper  layer.  Also,  0S  will  be 
greater  than  1 where  the  ceiling  surface  temperature  is  greater  than  TMAX.  This  is  possible,  for 
example,  during  times  of  reduced  fire  size  when  the  fire’s  near-ceiling  plume  temperature  is  reduced 
significantly,  perhaps  temporarily,  from  previous  values,  but  the  ceiling  surface,  heated  previously  to 
relatively  high  temperatures,  has  not  cooled  substantially.  Plots  of  0 per  Eq.  (32)  are  presented  in 
Figure  5 for  cases  when  ©s  is  < 0,  between  0 and  1,  and  > 0. 

In  order  to  calculate  TCJ(r,  z)  from  Eq.  (32),  it  is  necessary  to  determine  the  one  unknown 
parameter,  TMAX(r).  This  is  obtained  by  invoking  conservation  of  energy.  Thus,  at  an  arbitrary  r 
outside  the  stagnation  zone,  the  total  rate  of  radial  outflow  of  enthalpy  (relative  to  the  upper  layer 
temperature)  of  the  ceiling  jet  within  a small  subtended  angle  Ad  is  equal  to  the  fraction  A0/(27t)  of 
the  rate  of  enthalpy  flow  in  the  upper  layer  portion  of  the  plume,  Q\  less  the  integral  (from  the 
plume-ceiling  impingement  point  to  r)  of  the  flux  of  convective  heat  transfer  from  the  ceiling  jet  to 
the  ceiling  surface,  i.e., 
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When  0.2  < r/H: 


00 


2'7rr  /Pucp(tcj  ' Tu)VCJdz 
0 


- Q’  - 2rr  /q£EIL(r)rclr 
0 


(l-^CONv)Q’  (34) 


where  AEONv(r)  the  fraction  of  Q’  subtended  by  Ad  and  transferred  by  convection  to  the  ceiling 
from  the  point  of  ceiling  impingement  to  r,  i.e., 


^CONv(r)  ~ 


r 

[2^/qcEiLMrdr]/Q’ 


0 


(35) 


In  the  above,  Q’  has  been  calculated  previously  in  Eqs.  (11’).  Also,  for  particular  r’s  of  interest,  and 
with  q£EIL  evaluated  according  to  Eqs.  (13)-(19),  the  integral  on  the  right  hand  sides  of  Eqs.  (34)  and 
(35)  can  be  determined  numerically. 

The  integral  on  the  left  hand  side  of  Eq.  (34)  is  calculated  using  VCJ  of  Eqs.  (29)  and  (30)  and  TCJ 
of  Eqs.  (32)  and  (33).  From  this  the  desired  distribution  for  TMAX  is  found  finally  to  be  [13] 

(Tmax  - Tu)  = 2.6(1  - AioNv)(r/H)-0  8QftOTTu  - 0.090(TCEIL  - T^)  if  0.2  < r/H  (36) 

The  above  result  together  with  Eqs.  (32)  and  (33)  represent  the  desired  estimate  for  TCJ(r,z). 


Summary 

This  work  presented  a model  to  predict  the  instantaneous  rate  of  convective  heat  transfer  from  fire 
plume  gases  to  the  overhead  ceiling  surface  in  a room  of  fire  origin.  The  room  is  assumed  to  be  a 
rectangular  parallelopiped  and,  at  times  of  interest,  ceiling  temperatures  are  simulated  as  being 
uniform.  Also  presented  is  an  estimate  of  the  convective  heat  transfer,  due  to  ceiling-jet-driven  wall 
flows,  to  the  wall  surfaces.  The  effect  on  the  heat  transfer  of  the  location  of  the  fire  within  the  room 
is  taken  into  account.  Finally  presented  was  a model  of  the  velocity  and  temperature  distributions 
in  the  ceiling  jet. 

The  model  equations  were  used  to  develop  an  algorithm  and  associated  modular  computer  subroutine 
to  carry  out  the  indicated  heat  transfer  calculations.  The  subroutine  is  written  in  FORTRAN  77  and 
is  called  CEILHT.  A summary  flow  diagram  for  the  algorithm  is  presented  in  Figure  3.  The 
algorithm  and  subroutine  are  suitable  for  use  in  two-layer  zone-type  compartment  fire  model 
computer  codes.  CEILHT,  which  uses  the  two-dimensional  numerical  integration  subroutine  ADAPT 
[16]  is  presented  in  Appendix  B.  CEILHT  has  been  tested  for  a variety  of  instantaneous  fire 
environments  involving  a 107W  fire  in  a 8m  x 8m  x 4m  high  enclosure.  While  the  calculated  results 
were  plausible,  it  is  important  to  point  out  that  CEILHT  simulations  have  not  been  experimentally 
validated. 
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Nomenclature 


A\VALL,L,N  [A\VALL,U,n] 

area  of  lower-  [upper-]  layer  portion  of  wall  segment  N 

CP 

specific  heat  at  constant  pressure 

Cp 

9.115 

f(r/H) 

Eq.  (16) 

g 

acceleration  of  gravity 

H 

distance  of  continuation  source  below  ceiling,  Eq.  (10) 

B 

a characteristic  heat  transfer  coefficient,  Eq.  (8) 

hL 

heat  transfer  coefficient,  Eq.  (14) 

hST 

heat  transfer  coefficient  at  wall  stagnation  point,  Eq.  (21) 

M* 

fraction  of  Mplume  which  is  buoyant  in  upper  layer,  Eq.  (7) 

mplume 

mass  flow  rate  in  plume  at  ZLAYER  when  < ZLAYER 

Pr 

Prantl  number 

0’ 

strength  of  continuation  source,  Eq.  (IT) 

^CEIL 

rate  of  heat  transfer  to  ceiling,  Eq.  (2) 

^EQ 

dimensionless  strength  of  plume  at  ZLAYER 

0A 

dimensionless  strength  of  plume  at  ZCEIL 

^FCONV 

portion  of  fire  energy  release  rate  convected  in  plume 

QpLUME 

flow  of  enthalpy  in  plume  at  ZLAYER,  Eq.  (1) 

^ WALL 

rate  of  heat  transfer  to  all  wall  surfaces.  Eq.  (A-6) 

4ceil 

flux  of  rate  of  heat  transfer  to  ceiling,  Eq.  (13) 

4ceil,ave 

average  flux  of  rate  of  heat  transfer  to  ceiling,  Eq.  (2) 

^WALL 

flux  of  heat  transfer  to  the  wall 

9WALL,CRNN,L  M\VALL,CRNN,u] 

average  near  corner  N in  lower  [upper]  layer 
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4wall,l,ave  Wwall,u,ave] 

average  c^wall  Portion  of  wall  in  lower  [upper]  layer,  Eq.  (A-5) 

4wall,l,ave,n  Wwall,u,ave,n] 

average  to  Pordon  °f  wall  segment  N in  the  lower 

[upper]  layer 

3wall,norm,l  Wwall,norm,u] 

average  Qwall  to  l°wer"  [upper-]  layer  portion  of  wall  near 
vertical  line  through  ceiling-jet/wall  stagnation  point 

3\VALL,ST 

^wall  at  a ceiling-jet/wall  stagnation  point,  Eq.  (20) 

3wALL,STPN,L  WwALL,STPN,u] 

average  <)wall  to  portion  of  wall  in  the  lower-  [upper-]  layer 
along  the  line  passing  through  stagnation  point  N (STPN) 

ReH 

characteristic  Reynolds  number,  Eq.  (18) 

r 

distance  from  plume/ceiling  impingement  point 

rCRNN 

r of  corner  N (CRNN) 

rST 

r of  a stagnation  point 

rSTPN 

r of  stagnation  point  N (STPN) 

Tad 

adiabatic  ceiling  jet  temperature,  Eq.  (15) 

tceil 

temperature  of  ceiling 

TCj 

temperature  of  ceiling  jet 

Tl>  Ty 

temperature  of  lower  [upper]  layer 

twall 

average  temperature  of  the  wall  surfaces 

V 

characteristic  velocity  of  ceiling  jet,  Eq.  (30) 

Va 

velocity  of  ceiling  jet,  Eq.  (29) 

VMAX 

maximum  velocity  of  ceiling  jet,  Eq.  (30) 

X,Y,Z 

co-ordinates  in  room,  origin  at  a room  corner  on  the  floor 

XWALL’  YWALL>  and  ZCEIL 

X,  Y,  Z of  room  corner  diagonal  from  origin 

XFIRE’  Yfire,  and  ZFIRE 

X,  Y,  Z of  center  of  base  of  fire 

ZEQ 

Z of  equivalent  source  in  lower  layer,  Eq.  (6) 
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ZLAYER 

Z of  layer  interface 

ZSOURCE 

Z of  continuation  source,  Eq.  (12) 

Z 

distance  below  ceiling  surface,  Eq.  (28) 

a 

Tu/Tl 

5 

characteristic  ceiling  jet  thickness,  Eq.  (30)y 

0 

dimensionless  Tq,,  Eq.  (32) 

©s 

0 at  ceiling  surface 

^CONV 

Eq.  (35) 

vu 

kinematic  viscosity  of  air,  Eq.  (19) 

PL’  Pu 

average  density  of  lower  [upper]  layer 

a 

dimensionless  parameter,  Eq.  (9) 
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Appendix  A:  An  Algorithm  for  Estimating  the  Average  Rate  of  Ceiling-Jet-Driven  Convective  Heat 
Transfer  to  the  Upper  and  Lower  Layer  Walls 


1.  Use  Eqs.  (20)  and  (21)  to  calculate  <)wall  at  toP  the  waPs  near  the  points  STP1,  STP2, 
STP3,  and  STP4.  Then  use  Eqs.  (22)-(27)  to  calculate  the  average  ()wall  al°ng  the  upper-  and 
lower-layer  segments  of  the  vertical  lines  through  these  points. 

Refer  to  Figure  6.  Designate  the  radii  of  points  STP1,  STP2,  etc.  as  rSTP1,  rSTp->,  etc., 
respectively,  and  calculate  their  values: 


rSTPl  “ yfire 

rSTP2  = ^WALL  " XFIRE 

(A-l) 

rSTP3  “ YWALL  ’ YFIRE 
rSTP4  = XFIRE 

Use  rSTP1  as  rST  in  Eq.  (21)  and  as  r in  Eqs.  (15)  and  (17)  and  determine  the  corresponding 
value  of  4wall.  from  Eq.  (20).  Use  this  latter  value  as  wall  st  in  Eqs.  (22)-(27)  and  calculate 
the  desired  average  heat  transfer  along  the  upper  and  lower  line  segment  through  point  STP1. 
Designate  these  as  du/Air  stpi  ti  and  d™,  > ctpi  r - respectively.  Carry  out  analogous 
calculations  and  designates for the  STP2,  STP3,  and  STP4 

2.  Estimate  c^wall  at  the  top  °f  the  walls  near  the  four  corners  of  the  room  by  assuming  that  the 
effect  of  concentration  of  the  wall  flows  there  can  be  simulated  by  normal  ceiling-jet/wall 
by  using  Eqs.  (20)  and  (21).  For  example,  use  these  equations  with  rST  = 
1/2  to  estimate  the  heat  ransfer  near  the  top  of  the  walls  near  X = Y = 0. 
Then  use  Eqs.  (22)-(27)  to  calculate  the  average  <)wall  al°nS  the  upper-  and  lower-layer 
segments  of  the  vertical  lines  which  define  these  corners. 


impingement,  i.e., 
(xfire  + yfire) 


Refer  to  Figure  6.  Designate  the  radii  to  the  four  room  corners,  CRN1,  CRN2,  CRN3,  and 
CRN4  as  rCRN1,  rCRN2,  etc.,  respectively,  and  calculate  their  values: 


rCRNl  = (XFIRE  + YFIRe)1/2 

rCRN2  = t(XWALL  " XFIRe)~  + YFIRe]1/2 

rCRN3  = t(XWALL  ' XFIRe)~  + (YWALL  ' YFIRe)“]1/2 

rCRN4  = tXFIRE  + (YWALL  * YFIRe)“]1/2 


(A-2) 
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Replace  rST  in  Eq.  (21)  and  r in  Eqs.  (15)  and  (17)  by  rCRN1  and  determine  the  corresponding 
value  of  ()wall  fr°m  E<1-  (20).  Replace  <5wall,st  *n  ^qs.  (22)-(27)  by  this  latter  value  and 
calculate  the  desired  average  rates  of  heat  transfer  along  the  upper-  and  lower-layer  CRN1  line 
segment.  Designate  these  as  ()wall  crni  u anc*  3wall  crni  L>  respectively.  Carry  out 
analogous  calculations  and  designations  for  the  CRN2,  CRN3,  and  CRN4. 

3.  For  an  upper-  or  lower-layer  wall  segment  between  STP1,  STP2,  STP3,  or  STP4,  and  an 
adjacent  corner,  estimate  the  average  value  of  <5wall  as  average  °f  the  previously 
determined  average  values  along  the  vertical  bounding  "end-lines."  Calculate  this  average  for 
each  of  the  eight  upper-layer  segments  and  eight  lower-layer  segments. 

Refer  to  Figure  6.  Number  the  eight  pair  of  upper-  and  lower-layer  wall  segments  as  shown. 
Designate  the  average  rates  of  heat  transfer  to  upper-  and  lower-layer  segment  1 as 
4 wall  u ave  l anc*  ^wall  L ave  i>  respectively.  Desigate  analogous  rates  of  heat  transfer  to 
the  other  seven  pair  of  wall  segments.  Calculate: 


4wall,u,ave,i  = 
4wall,l,ave,i  = 

4 WALL, U,  AVE, 2 = 

4 WALL, L, AVE, 2 = 

4 WALL, U,  AVE, 3 = 

4 WALL, L,  AVE, 3 = 

^WALL,UAVE,4  = 
A WALL, L,  AVE, 4 = 

A\VALL,U,AVE,5  = 
^ WALL, L,  AVE, 5 = 

4wall,u,ave,6  = 

^ WALL, L, AVE, 6 = 

4waLL,U,AVE,7  = 
4wALL,L,AVE,7  = 
^ WALL, U,  AVE, 8 = 

3wALL,L,AVE,8  = 


MwaLL,CRN1,U 
(4  WALL, CRNI, L 
(AwaLL,CRN2,U 
(4waLL,CRN2,L 
(4wALL,CRN2,U 
(QwaLL,CRN2,L 
WwALL,CRN3,U 
(4\VALL,CRN3,L 
(^WALL,CRN3,U 
(^WALL,CRN3,L 
(^WALL,CRN4,U 
(3\VALL,CRN4,L 
WwALL,CRN4,U 
(^WALL,CRN4,L 
(4  WALL, CRNI, U 

(4wall,crni,l 


+ ^WALL,STPl,u)/2 
+ AwaLL,STP1,l)/2 
+ ‘iwALL.STPl.uV2 
+ 4wall,stpi,l)/2 
+ ^WALL,STP2,u)//2 
+ AwaLL,STP2,l)/2 
+ AwALL,STP2,u)y,2 
+ ^WALL,STP2,l)//2 
+ AwaLL,STP3,u)/2 
+ AwaLL,STP3,l)/2 
+ ^WALL.STPS.u)/2 

+ 4wall,stp3,l)/2 
+ Awall,stp4,u)/2 
+ 4wall,stp4,l)/2 

+ ^WALL,STP4,u)/2 
+ ^WALL,STP4,l)/2 


(A-3) 
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4.  For  each  of  the  upper-layer  segments  use  the  area  of  the  segment  and  the  previously  calculated 
average  cjwall  to  calculate  the  eight  contributions  to  the  total  rate  of  upper-layer  wall  heat 
transfer.  Sum  these  contributions  and  obtain  finally  <5wall  u ave’  the  average  flux  of  heat 
transfer  to  the  upper-layer  portions  of  the  walls.  Carry  out  analogous  calculations  for 
4 wall  l AVE’  the  average  flux  of  heat  transfer  to  the  lower-layer  portions  of  the  walls. 

Designate  the  area  of  upper-  and  lower-layer  wall  segment  N as  Awallun  and  Awallln, 
respectively,  and  calculate  these  values: 


AWALL,U,1 

= 

AWALL,U,6 

= 

awall,l,i 

= 

AWALL,L,6 

= 

AWALL,U,2 

= 

AWALL,U,5 

= 

AWALL,L,2 

= 

AWALL,L,5 

= 

AWALL,U,3 

= 

AWALL,U,8 

= 

AWALL,L,3 

= 

AWALL,L,8 

= 

AWALL,U,4 

= 

AWALL,U,7 

= 

AWALL,L,4 

= 

AWALL,L,7 

= 

Calculate 


xfire(zceil  ' zlayer) 

XFIREZLAYER 

(XWALL  * XFIRE)(ZCEIL  " ZLAYEr) 
(XWALL  " XFIRE)ZLAYER 

yfire(zceil  * zlayer) 
afirezlayer 

(ywall  • afire)(zceil  * zlayer) 

(^WALL  ’ AFIRe)ZLAYER 


o 

4wall,u,ave  = j|= (awall,u,n4  wall,  ave,  u,nV 


jf=|^WALL,U,N 


Awall,l,ave  “ 


f = [AWAll,L,N^  WALL, AVE, L,  nV 


8 

fl=j^WALL,L,N 


(A-4) 


(A-5) 


5.  Add  the  rates  of  heat  transfer  to  all  sixteen  upper  and  lower  wall  segments  and  obtain,  OwalL’ 
the  total  rate  of  heat  transfer  to  the  wall  surfaces. 


0 


WALL  ~ 3wALL,U,AVE 


8 „ Y 

§=1AWALL,U,N  + 4 WALL, L,  AVE  AWALL,L,N 


(A-6) 
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Appendix  B:  The  Subroutine  CEELHT  for  Calculating  the  Convective  Heat  Transfer  in  a Two-Layer 
Fire  Environment  From  Fire  Plume  Gases  to  a Uniform  Temperature  Ceiling  and  to 
Uniform  Temperature  Walls. 


Input  and  Output  Variables 

The  CEILHT  input  and  output  variables  with  definitions  and  units  are  listed  below.  Except  for  EPS, 
these  are  taken  from  the  descriptions  presented  in  the  body  of  the  text  and  in  Appendix  A 

OUTPUT 

Qceil 

rate  of  heat  transfer  to  ceiling  [W] 

^WALL 

rate  of  heat  transfer  to  all  wall  surfaces  [W] 


3ceil,ave 

average  flux  of  rate  of  heat  transfer  to  ceiling  [W/m“] 
4wall,l,ave  Wwall,u,ave] 

average  heat  flux  to  portion  of  wall  in  lower  [upper]  layer  [W/m2] 


INPUT 

MpLUME 


mass  flow  rate  in  plume  at  2^^  ayfr  tf  ^ ^layer  [kg/s] 


0 


FCONV 


portion  of  fire  energy  release  rate  convected  in  plume  [W] 
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tceil 


average  temperature  of  ceiling  [K] 


TL’  Ty 

temperature  of  lower  [upper]  layer  [K] 

twall 

average  temperature  of  the  wall  surfaces  [K] 

XwalL’  YwalL’  and  ZCEIL 

X,  Y,  Z of  room  corner  diagonal  from  origin  [m] 

XFIRE’  YFIRE’  and  ZFIRE 

X,  Y,  Z of  center  of  base  of  fire  [m] 


Z LAYER 


Z of  layer  interface  [m] 


PL’  Pu 

density  of  lower  [upper]  layer  [kg/m3] 

Cross-Reference  Between  Nomenclature  and  Subroutine  Variables 

Below  is  a cross-reference  between  variables  defined  in  the  Nomenclature  section  and  variables  used 
in  the  subroutine. 

NOMENCLATURE  VARIABLES  SUBROUTINE  VARIABLES 

AWall,l,n>  Awall,U)N  AWU(N),  AWU(N)  [m2] 

Cp  CP  (use  1000.)  [W-s/(kg-K)] 
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CT 

f(r/H) 

g 

H 

fi 

hL 

hST 

M* 

mplume 

Pr 

0’ 

OcEIL 

Oeq 

Oh 

Ofconv 

OpLUME 

OwALL 

Aceil 

Aceil,ave 

A wall 

A\VALL,CRNN,L’  Awall,crnn,u 
AwALL,L,AVE’  Awall,u,ave 
AwaLL,L,AVE,N’  Awall,u,ave,n 
AwaLL,NORM,L’  Awall,norm,u 


CT  (use  9.115) 

F 

G (use  9.8)  [m/s2] 

H [m] 

HTCT  [W/(m2-K)] 

HTCL  [W/(m2-K)] 

HTCST  [W/(m2-K)j 
MFRAC 
MPLUME  [m/s] 

PR 

QCONT  [W] 

QCEIL  [W] 

QEQ  [W] 

QH 

QCONV  [W] 

QPLUME  [W] 

QWALL  [W] 

QFCLG  [W/m2] 

QFCLGA  [W/m2] 

QFW  [W/m2] 

QFWCL(N),  QFWCU(N)  [W/m2] 
QFWLA,  QFWUA  [W/m2] 
QFWLAN(N),  QFWUAN(N)  [W/m2] 
QFWNL,  QFWNU  [W/m2] 
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4wall,st 

QFWST  [W/m2] 

4waLL,STPN,L’  4wall,stpn,u 

QFWSL(N),  QFWSU(N)  [W/m2 

ReH 

RE 

r 

R [m] 

rCRNN 

RC(N)  [m] 

rST 

RST  [m] 

rSTPN 

RS(N)  [m] 

Tad 

TAD  [K] 

tceil 

TC  [K] 

Tcj 

TCJ  [K] 

tl>  Tu 

TL,  TU  [K] 

twall 

TW  [K] 

V 

VCH  [m/s] 

Vcj 

VCJ  [m/s] 

VMAX 

VMAX  [m/s] 

X,Y,Z 

X,  Y,  Z [m] 

XWALL’  Ywall,  and  ZCEIL 

XW,  YW,  ZC  [m] 

^FIRE’  T"nRE,  and  ZnRE 

XF,  YF,  ZF  [m] 

ZEQ 

ZEQ  [m] 

ZLAYER 

ZLAY  [m] 

ZSOURCE 

ZS  [m] 

Z 

z [m] 

a 

ALPHA 

8 

DEL  [m] 
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0 


THETA 


©s 

THETAS 

^•CONV 

LAMC 

VU 

ANU  [m2/s] 

PL’  Pu 

RHOL,  RHOU  [kg/m3] 

a 

SIGMA 

The  Subroutine  CEILHT 

On  the  following  pages  is  a listing  of  the  subroutine  CEILHT. 
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SUBROUTINE  CEILHT( 

I MPLUME , QCONV , ATC , TL , TU , TW , XW , YW , ZC , AXF , 

I AYF , ZF , ZLAY , RHOL , RHOU , EPS , 

0 QCEIL , QFCLGA , QWALL , QFWLA , QFWUA) 

IMPLICIT  DOUBLE  PRECISION  (A-H.M.O-Z) 

INTEGER  RULCLS , MAXPTS , MINPTS 


C*** 

c*** 


Q ~k  -k  ~k  'k  ~k  ~k  -k  k k k k k 'A'  k k k -A' k k k '.V  k k k k k k k k 'V  k k '-V  k k k k k k k k k k k k k k k k -/V  k k - V k k k k k k k k k k k k Vc  k k 

C*******WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING 
C***  This  subroutine  and  the  function  QFCLG  used  in  its  called 
subroutines  use  the  common  blocks  AINTCH 

***WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING 
C***  The  algorithm  on  which  this  subroutine  is  based  has  not  been 
C***  validated  experimentally 

Q "■rC  /C  VC"  VC"  VC"  VC  VC  VC  VC  /f  "X"  "a  /\  "/C  VC  /v  /i"  "a  'C  /\  /wC  VC  X VC  /C  VC  a"'1*  VC  X a"  VC  VC  1C  /v  VC  a"  VC  'f /C  VC  Iwf  /\  VC"  "/C/i  'C  "id  1C  1C  /C  VC  1C  /v 

THIS  SUBROUTINE  CALCULATES  CONVECTIVE  HEAT  TRANSFER  TO 
THE  UNIFORM  TEMPERATURE  CEILING  ABOVE  A FIRE  IN  A PARALLEL- 
OPIPED  ROOM  WITH  A TWO - LAYER  FIRE  ENVIRONMENT.  ALSO  CALCU- 
LATED IS  THE  TOTAL  RATE  OF  HEAT  TRANSFER  TO  THE  WALLS  AND 
THE  AVERAGE  FLUX  TO  THE  UPPER  AND  LOWER  PORTIONS  OF  THE 


WALLS 


INPUT 

c*** 

c 

Qkkk 

c*** 

EPS 

c*** 

c*** 

MPLUME 

c*** 

Qkkk 

QCONV 

c*** 

TC 

c*** 

TL,  TU 

c*** 

TW 

c*** 

C'k'k'k 

Q XXX 

XW , YW , ZC 

C 'k  i’c  "k 

XF , YF , ZF 

Qkkk 

ZLAY 

c*** 

c 

RHOL, RHOU 

C'k'-k'k 

OUTPUT 

c*** 

c 

c*** 

QCEIL 

c*** 

Qkkk 

QFCLGA 

Qkkk 

QWALL 

c*** 

C 'k  "k  ~k 

C'-k'k'k 

QFWLA , QFWUA 

RELATIVE  ACCURACY  FOR  COMPUTING  THE  INTEGRAL 
OF  THE  CEILING  HEAT  TRANSFER  FLUX 

MASS  FLOW  RATE  IN  THE  PLUME  AT  ZLAY  IF  ZF  < 
ZLAY  [KG/S] 

PORTION  OF  FIRE  ENERGY  RELEASE  RATE  CONVECT- 
ED  IN  PLUME  [W] 

AVERAGE  TEMPERATURE  OF  CEILING  [K] 

TEMPERATURE  OF  LOWER,  UPPER  LAYER  [K] 

AVERAGE  TEMPERATURE  OF  WALL  SURFACES  [K] 
CO-ORDINATES  OF  ROOM  CORNER  DIAGONALLY  OPPO- 
SITE TO  ORIGIN  OF  AXES  (ORIGIN  IS  AT  A CORN- 
ER ON  THE  FLOOR,  WITH  X,  Y AZES  ALONG  WALL/ 

FLOOR  JUNCTION  AND  Z AXES  UPWARD)  [M] 

CO-ORDINATES  OF  CENTER  OF  BASE  OF  FIRE  [M] 
ELEVATION  ABOVE  FLOOR  OF  LAYER  INTERFACE  [M] 
DENSITY  OF  LOWER,  UPPER  LAYER  [KG/M**3] 


RATE  OF  HEAT  TRANSFER  TO  CEILING  [W] 

AVERAGE  FLUX  OF  HEAT  TRANSFER  TO  CEILING 
[W/M**2 ] 

RATE  OF  HEAT  TRANSFER  TO  WALL  [W] 

AVERAGE  FLUX  OF  HEAT  TRANSFER  TO  LOWER  AND 
UPPER  PORTIONS  OF  THE  WALL  SUFACES  [W/M**2] 
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SOME  OTHER  DEFINITIONS  AND  FIXED  INPUT  IN  THIS  SUBROUTINE 


c*** 

c*** 

C.kkk 

SOME  OTHE 

ALPHA 

Qkkk 

AWL(N) 

Q kkk 

AWU(N) 

Qkkk 

CP 

Qkkk 

Qkkk 

CT 

Qkkk 

G 

Q kkk 

H 

Qkkk 

Qkkk 

Qkkk 

Qkkk 

QCONT 

Qkkk 

QEQ 

C kkk 

QH 

Qkkk 

QFWCL(N) 

Qkkk 

Qkkk 

QFWCU(N) 

Qkkk 

Qkkk 

QFWLAN (N) 

Qkkk 

Qkkk 

QFWUAN (N) 

Qkkk 

Qkkk 

QFWSL(N) 

Qkkk 

Qkkk 

Qkkk 

QFWSU(N) 

Qkkk 

Qkkk 

Qkkk 

QFWST 

C kkk 
Qkkk 

PR 

Qkkk 

RC  (N) 

Qkkk 

Qkkk 

RS  (N) 

Q'k'k'k 

THT 

Qkkk 

Qkkk 

ZEQ 

Qkkk 

Qkkk 

ZS 

Qkkk 

TU/TL 

AREA  OF  LOWER - LAYER  PORTION  OF  WALL  SEGMENT  N 
AREA  OF  UPPER- LAYER  PORTION  OF  WALL  SEGMENT  N 
SPECIFIC  HEAT  OF  AIR  AT  CONSTANT  PRESSURE 
[KJ/(KG*K) ] 

9.115,  CONSTANT  IN  POINT  SOURCE  PLUME  EQN . 

9.8,  ACCELERATION  OF  GRAVITY  [M/S**2] 

IF  LAYER  INTERFACE  IS  BELOW  CEILING:  DISTANCE 
OF  CONTINUATION  SOURCE  BELOW  CEILING;  IF  LAY- 
ER INTERFACE  IS  AT  CEILING:  DISTANCE  OF  EQUI- 
VALENT SOURCE  BELOW  CEILING  [M] 

STRENGTH  OF  CONTINUATION  SOURCE  [W] 

DIMENSIONLESS  STRENGTH  OF  PLUME  AT  ZLAY 
DIMENSIONLESS  STRENGTH  OF  PLUME  AT  ZC 
AVERAGE  HEAT  FLUX  TO  WALL  NEAR  CORNER  N IN 
THE  LOWER  LAYER  [W/M**2] 

AVERAGE  HEAT  FLUX  TO  WALL  NEAR  CORNER  N IN 
THE  UPPER  LAYER  [W/M**2] 

AVERAGE  HEAT  FLUX  TO  PORTION  OF  WALL  SEGMENT 
N IN  THE  LOWER  LAYER  [W/M**2] 

AVERAGE  HEAT  FLUX  TO  PORTION  OF  WALL  SEGMENT 
N IN  THE  UPPER  LAYER  [W/M**2] 

AVERAGE  HEAT  FLUX  TO  PORTION  OF  WALL  IN  THE 
LOWER  LAYER  ALONG  THE  LINE  PASSING  THROUGH 
STAGNATION  POINT  N [W/M**2] 

AVERAGE  HEAT  FLUX  TO  PORTION  OF  WALL  IN  THE 
UPPER  LAYER  ALONG  THE  LINE  PASSING  THROUGH 
STAGNATION  POINT  N [W/M**2] 

HEAT  TRANSFER  FLUX  TO  WALL  AT  A WALL/CEILING 
-JET  STAGNATION  POINT  [W/M**2] 

PRANDTL  NUMBER 

DISTANCE  FROM  PLUME/CEILING  IMPINGEMENT  POINT 
TO  CORNER  N [M] 

DISTANCE  FROM  PLUME/CEILING  IMPINGEMENT  POINT 
TO  WALL  STAGNATION  POINT  N [M] 

TEMPERATURE  OF  AMBIENT  IN  UNCONFINED  CEILING 
HEAT  TRANSFER  PROBLEM,  EITHER  TL  OR  TU  [K] 

ELEVATION  ABOVE  FLOOR  OF  EQUIVALENT  SOURCE 
IN  LOWER  LAYER  [M] 

ELEVATION  ABOVE  FLOOR  OF  CONTINUATION  SOURCE 
[M] 


DIMENSION  A( 2 ) , B ( 2 ) 

PARAMETER  (LENWRK=10000) 

DIMENSION  WRKSTR ( LENWRK) 

EXTERNAL  QFCLG 

PARAMETER (PI-3 . 141592D0 , G=9 . 8D0 , CT=9 . 115D0 , CP=1000 . DO , 
+ PR=0 . 7D0 ) 

DIMENSION  AWL ( 8 ) ,AWU(8) ,QFWCL(4) ,QFWCU(4) , QFWLAN ( 8 ) , 

+ QFWUAN ( 8 ) , QFWSL(4) ,QFWSU(4) ,RC(4) ,RS(4) 
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c 


c 


c 


c 


COMMON/AINTCH/H , HTCT , THT , THTQHP , Cl , C2 , C3 , XF , YF , TC 

XF=AXF 

YF=AYF 

TC=ATC 

IF(ZF. LT . ZLAY)THEN 

fire  is  below  layer  interface************************ 
QEQ=(0.21D0*QCONV/(CP*TL*MPLUME))**(3.D0/2.D0) 
ZEQ=ZLAY- (QCONV/(QEQ*RHOL*CP*TL*DSQRT(G) ) )** 

+ (2. DO/5. DO) 

I F ( ZLAY . LT . ZC ) THEN 

layer  interface  is  below  ceil ing************* 
ALPHA=TU/TL 
IF  ( ALPHA. LE. 1. DO) THEN 
MFRAC=1 . DO 
QC0NT=QC0NV 
ZS=ZEQ 
THT=TU 
RHOHT=RHOU 


+ 

+ 


+ 


+ 

+ 

+ 

+ 


ELSE 

S I GMA= ( 1 . DO  - ALPHA+CT* 

(QEQ**(2 . DO/3 . DO) )) /(ALPHA- 1 . DO) 
IF ( SIGMA. GT . 0 . DO) THEN 
SSQ=SIGMA**2.D0 
T0P=1 . 04599D0*SIGMA+ 

0. 360391D0*SSQ 

B0TT0M=1 . DO+1 . 37748D0*SIGMA+ 

0.360391D0*SSQ 

MFRAC=TOP/BOTTOM 

QCONT-QCONV*SIGMA*MFRAC/ 

( 1 . DO+SIGMA) 
ZS=ZLAY - (ZLAY-ZEQ)* 

(ALPHA** ( 3 . DO/5 . DO ) ) * 
(MFRAC** ( 2 . DO/ 5 . DO) )* 

( ( (1. DO+SIGMA) /SIGMA) 
**(1 . DO/5 . DO) ) 

THT=TU 

RHOHT=RHOU 


ELSE 

no  plume  penetration********** 

MFRAC=0 . DO 

QCEIL=0 . DO 

QFCLGA=0 . DO 

QWALL=0 . DO 

QFWLA=0 . DO 

QFWUA=0 . DO 

RETURN 

ENDIF 


END  IF 

ELSE 

layer  interface  is  at  ceiling***************** 

QCONT=QCONV 

ZS=ZEQ 
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THT=TL 

RHOHT=RHOL 


ENDIF 


ELSE 


fire  is  at  or  above  layer  interface*1* 
IF( ZF . LT . ZC) THEN 

fire  is  below  ceiling*********1* 

QCONT=QCONV 

ZS=ZF 

THT=TU 

RHOHT=RHOU 

ELSE 


fire  is  at  ceiling************ 


Vc~A-b'wV-, 


QCEIL=0 . DO 
QFCLGA=0 . DO 
QWALL=0 . DO 
QFWLA=0 . DO 
QFUUA=0 . DO 
RETURN 

ENDIF 

ENDIF 

H=ZC-ZS 

SQRTGH= ( G*H ) ** ( 1 . DO/2 . DO ) 
QH=QCONT/(RHOHT*CP*THT*SQRTGH* (H**2 . DO) ) 

QHP= (QH** ( 1 . DO/3 . DO ) ) 

HTCT=RHOHT*CP*SQRTGH*QHP 

ANU=(0 . 04128D-7*(THT**(5 . DO/2 . DO) ) ) /(THT+110 . 4D0) 

RE=SQRTGH*H*QHP/ANU 

THTQHP-THT* (QHP**2 . DO ) 

PRP=PR** ( 2 . DO/3 . DO) 

Cl=8 . 82D0/(DSQRT (RE)*PRP) 

C2=5 . DO  - 0 . 2 84D0* (RE**0 . 2D0) 

C3=0 . 283032655 D0/( (RE**0 . 3D0)*PRP) 

C4=0 . 94D0/( (RE**0 . 42D0)*PR) 

NDIM  = 2 
MAXPTS  = 25000 

RULCLS  = 2 **ND IM+2  *ND IM**  2+ 6 *ND IM+ 1 

ITEMP1  = ( 2*NDIM+3 ) * ( 1+MAXPTS /RULCLS ) /2 
ITEMP  = 2 **ND IM+  2*ND IM**  2+6  *ND IM+ 1 
QCEIL=0 . 

DO  5 1=1,4 
MINPTS  = 10 
IF(I.EQ. 1)THEN 


A(  1) 
A(2). 
B(l) 
B ( 2 ) 


= 0. 
= 0. 
= XF 
= YF 


ELSE 


IF(I . EQ. 2) THEN 
A( 1)  = XF 
A(2)  = 0. 
B(l)  = XW 
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ELSE 


B(2)  = YF 


IF( I . 

EQ . 3 ) 

THEN 

A(  1) 

= XF 

A(2) 

= YF 

B ( 1 ) 

= XW 

B(2) 

= YW 

ELSE 

A(l) 

= 0. 

A(2) 

= YF 

B(l) 

= XF 

B (2 ) 

= YW 

ENDIF 

END  IF 

ENDIF 

CALL  ADAPT (NDIM , A, B , MINPTS , MAXPTS , QFCLG , EPS , RELERR , LENWRK , 
* WRKSTR, ADDQCEIL, IFAIL) 

QCEIL=QCEIL+ADDQCEIL 
5 CONTINUE 

IF(IFAIL.GT.O)  QCEIL=0 . DO 
QFCLGA=QCEIL/(XW*YW) 

C***  Now  calculate  wall  heat  transfer: 

C***  Step  1. 

C***  Calculate  radii  at  wall  stagnation  points: 

RS ( 1 ) =YF 
RS (2)=XW-XF 
RS(3)=YW-YF 
RS (4)=XF 

RMIN=MIN(RS ( 1) , RS ( 2 ) , RS ( 3 ) , RS ( 4 ) ) 

IF  (RMIN.LT.0.2*(ZC-ZF))THEN 

WRITE ( 6 , *) ' IMPINGEMENT  POINT  IS  CLOSER  TO  A WALL  THAN  0.2  OF 
+ ' FIRE-TO-CEILING  DISTANCE' 

ENDIF 

C***  Calculate  average  heat  transfer  fluxes  to  lower  and  upper 
C***  walls  along  the  vertical  lines  passing  through  the  four  wall/ 
C***  ceiling- jet  stagnation  points: 

DO  10  N=1 , 4 
RDH=RS (N)/H 

CALL  SQFWST(RDH , H , C4 , THT , HTCT , THTQHP , TW , QFWSL(N) ,QFWSU(N) , 

+ ZC , ZLAY) 

10  CONTINUE 
C***  Step  2. 

C***  Calculate  radii  at  room  corners: 

RC ( 1 ) =DSQRT (XF**2 . D0+YF**2 . DO ) 

RC ( 2 )=DSQRT ( (XW-XF)**2 . D0+YF**2 . DO) 

RC ( 3 ) =DSQRT ( (XW-XF)**2 . D0+(YW-YF)**2 . DO ) 

RC(4)=DSQRT(XF**2 . D0+(YW-YF)**2 . DO) 

C***  Calculate  average  heat  transfer  fluxes  to  lower  and  upper 
C***  walls  along  the  vertical  lines  passing  through  the  four  room 
C***  corners  by  assuming  that  the  heat  transfer  there  is  as  along 
C***  a line  passing  through  a point  of  normal  ceiling- j et/wall  im- 
C***  pingement. 
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DO  20  N=1 , 4 
RDH=RC(N)/H 

CALL  SQFWST (RDH , H , C4 , THT , HTCT , THTQHP , TW , QFWCL(N) , QFWCU (N) , 

+ ZC , ZLAY) 

20  CONTINUE 

C***  Step  3. 

C***  Calculate  the  average  heat  transfer  fluxes  to  the  lower  and 
C***  upper  portions  of  the  eight  wall  segments  bounded  by  the  room 
C***  corners  and  the  the  vertical  lines  passing  through  the  points 
C***  of  normal  wall/ceiling- j et  impingement. 
QFWUAN(1)=(QFWCU(1)+QFWSU(1) )/2 . DO 
QFWIAN ( 1 ) = ( QFWCL ( 1 ) +QFWS L ( 1 ) ) /2 . DO 
QFWUAN(2)=(QFWCU(2)+QFWSU(1) )/2 . DO 
QFWLAN ( 2)=(QFWCL( 2 )+QFWSL( 1) ) /2 . DO 
QFWUAN ( 3 )=(QFWCU ( 2)+QFWSU ( 2) )/2 . DO 
QFWLAN ( 3 )=(QFWCL( 2 )+QFWSL( 2 ) )/2 . DO 
QFWUAN (4) =(QFWCU ( 3 )+QFWSU ( 2 ) )/2 . DO 
QFWLAN(4)=(QFWCL(3)+QFWSL(2) )/2 . DO 
QFWUAN ( 5 ) = ( QFWCU ( 3 ) +QFWSU ( 3 ) ) /2 . DO 
QFWLAN(5)=(QFWCL(3)+QFWSL(3) )/2 .DO 
QFWUAN ( 6 ) = ( QFWCU ( 4 ) +QFWSU ( 3 ) ) /2 . DO 
QFW1AN ( 6 ) = ( QFWCL ( 4 ) +QFWS L ( 3 ) ) / 2 . DO 
QFWUAN ( 7 )=(QFWCU(4)+QFWSU(4) ) /2 . DO 
QFWLAN ( 7 )=(QFWCL(4)+QFWSL(4) ) /2 . DO 
QFWUAN ( 8 ) = ( QFWCU ( 1 ) +QFWSU ( 4 ) ) /2 . DO 
QFWLAN(8)=(QFWCL(1)+QFWSL(4) )/2 . DO 
C***  Step  4. 

C***  For  each  of  the  upper  layer  segments  use  the  area  of  the  seg- 
C***  ment  and  the  previously  calculated  average  heat  transfer 
C***  flux  to  calculate  the  eight  contributions  to  theto  the  total 
C***  rate  of  upper- layer  wall  heat  transfer.  Sum  these  contribu- 
C***  tions  and  obtain  finally  the  average  rate  of  heat  transfer  to 
C***  the  upper- layer  portions  of  the  walls.  Carry  out  analogous 
C***  calculations  for  the  lower  wall  surfaces.  Add  rates  of  heat 
C***  transfer  to  all  16  wall  surface  segments  and  obtain  total 
C***  rate  of  heat  transfer  to  the  wall. 

AWL ( 1 ) =XF*ZLAY 
AWU ( 1 ) -XF* ( ZC - ZLAY) 

AWL ( 2 ) = (XW - XF ) -ZLAY 
AWU ( 2 ) = (XW - XF ) * ( ZC - ZLAY ) 

AWL(3)=YF*ZLAY 
AWU ( 3 ) =YF* ( ZC - ZLAY) 

AWL ( 4 ) = ( YW - YF ) *ZLAY 
AWU(4)=(YW-YF)*(ZC-ZLAY) 

AWL ( 5 ) =AWL ( 2 ) 

AWU ( 5 ) =AWU ( 2 ) 

AWL ( 6 ) =AWL ( 1 ) 

AWU ( 6 ) =AWU ( 1 ) 

AWL ( 7 ) =AWL ( 4 ) 

AWU ( 7 ) =AWU ( 4 ) 

AWL ( 8 ) =AWL ( 3 ) 


33 


AWU  ( 8 ) =AWU  ( 3 ) 

SUMAQL=0 . DO 
SUMAQU=0 . DO 
SUMAL=0 . DO 
SUMAU=0 . DO 
DO  30  N=1 , 8 

SUMAQL=AWL ( N ) *QFWLAN ( N ) +SUMAQL 
SUMAQU=AWU ( N ) *QFWUAN ( N ) +SUMAQU 
SUMAL=AWL ( N ) +SUMAL 
SUMAU=AWU ( N ) +SUMAU 
30  CONTINUE 

IF ( SUMAL . LE „ 0 . DO ) THEN 
QFWLA=0 . DO 

ELSE 

QFWLA=SUMAQL/SUMAL 

ENDIF 

I F ( SUMAU . LE , 0 . DO ) THEN 
QFWUA=0 „ DO 

ELSE 

Q FWUA= SUMAQU/ S UMAU 

ENDIF 

QWALL= S UMAQL+  SUMAQU 

RETURN 

END 

C ~k  'k  k ~k  'k  'k  k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k ~'c  k k k k k -A' 


Qkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk  kkkkkkkk k kkkk  k kkkkkkkkkkkkkk  kkk 

SUBROUTINE  SQFWST (RDH , H , C4 , THT , HTCT , THTQHP , TW , QFWLOW , QFWUP , 

+ ZC , ZLAY) 

C***  Calculate  average  heat  transfer  fluxes  to  lower  and  upper 
Q'k-k'k  walls  along  a vertical  line  passing  through  a wall/ceiling - 
C***  jet  stagnation  point 

IMPLICIT  DOUBLE  PRECISION  (A-H.M.O-Z) 
F-(1.D0-1.1D0*(RDH**0.8D0)+0.808D0*(RDH**1.6D0))/ 

+ ( 1 . DO  - 1 . IDO* (RDH**0 . 8D0 ) +2 . 2D0*(RDH**1 . 6D0 ) + 

+ 0 . 6 9D0* ( RDH**2 . 4D0 ) ) 

IF(RDH . LT . 0 . 2D0)THEN 

TADDIM=10 . 22D0 - 14 . 9D0*RDH 

ELSE 

TADDIM=8.39D0*F 

ENDIF 

HTCL-C4*HTCT/RDH 
TAD=TADDIM*THTQHP+THT 
QFWST-HTCL* ( TAD - TW) 

IF(ZC.LE.0.8D0*H)THEN 

QFWLOW= (QFWST/2 . DO ) * ( 1 . 6D0*H - ( ZC - ZLAY) - ZC ) / ( 0 . 8D0*H) 
QFWUP=QFWST* ( 1 . DO - ( ZC - ZLAY) / ( 1 . 6D0*H) ) 

ELSE 

IF( (ZC-ZLAY) .GT.O. 8D0*H)THEN 
QFWLOW=0 . DO 

QFWUP= (QFWST/2 . DO ) * ( ZC - 0 . 8D0*H) / ( ZC - ZLAY) 

ELSE 


34 


QFWLOW= (QFWST/2 . DO ) * ( 0 . 8DO*H - ( ZC - ZLAY) ) /ZLAY 
QFWUP=QFWST*(1.D0- (ZC-ZLAY)/(1.6D0*H) ) 

END  IF 


ENDIF 

RETURN 

END 

Q ~k  -k  ~k  ~k  -k  ~k  ~k  ~k  -k  -k  ~k  -k  ~.V  it  -k  -k  'k  ~k  -k  k ~k  ~k  k k k k k k k k k k k k it  k k k k k k k k k k k k k -,V  i'<  -A'  Vc  -A"  - V Vc  Vc  -A-  -A-  Vc  Vc  v’c  Vc  -A-  Vc  V.-  Vc  V' 
Vf  'A -V  -A  A' A'  A-  A A A A A A A A A A A A A A A A A A A A A A A A A-  A A A'  A A A A A A-  A A A A A A A A A A A A A A A A A A A A A A 
CAAAAAAAAAAAAAAAAAAAAAA'AAAAAA'AAAAAAAAAAAAAAA'A'AAA'AA  aaa-aaaaaaa-aaaaaaaaa 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C 

C 

C 

C 

C 

C 


SUBROUTINE  ADAPT (NDIM , A , B , MINPTS , MAXPTS , FUNCTN , EPS , RELERR , LENWRK , 
* WRKSTR, FINEST, IFAIL) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

AAABEGIN  PROLOGUE  ADAPT 

ADAPTIVE  MULTIDIMENSIONAL  INTEGRATION  SUBROUTINE 

AUTHOR:  A.  C.  GENZ,  Washington  State  University- 
19  March  1984 

AAAAAAAAAAAAAA  PARAMETERS  FOR  ADAPT  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

aaaaa  INPUT  parameters 

NDIM  NUMBER  OF  VARIABLES,  MUST  EXCEED  1,  BUT  NOT  EXCEED  20 

A REAL  ARRAY  OF  LOWER  LIMITS,  WITH  DIMENSION  NDIM 

B REAL  ARRAY  OF  UPPER  LIMITS,  WITH  DIMENSION  NDIM 

MINPTS  MINIMUM  NUMBER  OF  FUNCTION  EVALUATIONS  TO  BE  ALLOWED. 

ON  THE  FIRST  CALL  TO  ADAPT  MINPTS  SHOULD  BE  SET  TO  A 
NON  NEGATIVE  VALUE.  (CAUTION...  MINPTS  IS  ALTERED  BY  ADAPT) 
IT  IS  POSSIBLE  TO  CONTINUE  A CALCULATION  TO  GREATER  ACCURACY 
BY  CALLING  ADAPT  AGAIN  BY  DECREASING  EPS  (DESCRIBED  BELOW) 
AND  RESETTING  MINPTS  TO  ANY  NEGATIVE  VALUE. 

MINPTS  MUST  NOT  EXCEED  MAXPTS. 

MAXPTS  MAXIMUM  NUMBER  OF  FUNCTION  EVALUATIONS  TO  BE  ALLOWED, 

WHICH  MUST  BE  AT  LEAST  RULCLS , WHERE 
RULCLS  = 2**NDIM+2*NDIM**2+6*NDIM+1 


FOR  NDIM  =23456789  10 

MAXPTS  >=  25  45  73  113  173  269  433  729  1285 

A suggested  value  for  MAXPTS  is  100  times  the  above  values. 


FUNCTN  EXTERNALLY  DECLARED  USER  DEFINED  FUNCTION  TO  BE  INTEGRATED. 
IT  MUST  HAVE  PARAMETERS  (NDIM.Z),  WHERE  Z IS  A REAL  ARRAY 
OF  DIMENSION  NDIM. 

EPS  REQUIRED  RELATIVE  ACCURACY 

LENWRK  LENGTH  OF  ARRAY  WRKSTR  OF  WORKING  STORAGE,  THE  ROUTINE 
NEEDS  ( 2*NDIM+3 ) * ( 1+MAXPTS /RULCLS ) /2  FOR  LENWRK  IF 
MAXPTS  FUNCTION  CALLS  ARE  USED. 

FOR  GUIDANCE,  IF  YOU  SET  MAXPTS  TO  100*RULCLS  (SEE  TABLE 
ABOVE)  THEN  ACCEPTABLE  VALUES  FOR  LENWRK  ARE 

FOR  NDIM  = 2345678  9 

LENWRK  = 357  561  1785  3417  6681  13209  26265  52377 


*****  OUTPUT  PARAMETERS 

MINPTS  ACTUAL  NUMBER  OF  FUNCTION  EVALUATIONS  USED  BY  ADAPT 
WRKSTR  REAL  ARRAY  OF  WORKING  STORAGE  OF  DIMENSION  (LENWRK) . 
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C RELERR  ESTIMATED  RELATIVE  ACCURACY  OF  FINEST 
C FINEST  ESTIMATED  VALUE  OF  INTEGRAL 

C I FAIL  IFAIL=0  FOR  NORMAL  EXIT,  WHEN  ESTIMATED  RELATIVE  ACCURACY 

C RELERR  IS  LESS  THAN  EPS  WITH  MAXPTS  OR  LESS  FUNCTION 

C CALLS  MADE . 

C IFAIL=1  IF  MAXPTS  WAS  TOO  SMALL  FOR  ADAPT  TO  OBTAIN  THE 

C REQUIRED  RELATIVE  ACCURACY  EPS.  IN  THIS  CASE  ADAPT 

C RETURNS  A VALUE  OF  FINEST  WITH  ESTIMATED  RELATIVE 

C ACCURACY  RELERR. 

C IFAIL=2  IF  LENWRK  TOO  SMALL  FOR  MAXPTS  FUNCTION  CALLS.  IN 

C THIS  CASE  ADAPT  RETURNS  A VALUE  OF  FINEST  WITH 

C ESTIMATED  ACCURACY  RELERR  USING  THE  WORKING  STORAGE 

C AVAILABLE,  BUT  RELERR  WILL  BE  GREATER  THAN  EPS. 

C IFAIL=3  IF  NDIM  ) 2,  NDIM  \ 20,  MINPTS  \ MAXPTS, 

C OR  MAXPTS  ) RULCLS . 

Q -k  Vc  -k  k k k kkk  k k k k k k '.V  k k k k k k k k k k k k k k k k k k k k k -k  ~k  ~k  -k  -!<  -k  ~k  'k  -k  ~k  ~k  -k  ~k  ~k -k  -k  -A'  -A'  -A-  -A-  -,V  -k  -k  -k  -k  'k  -,V  Vc 

C'A-*-*END  PROLOGUE  ADAPT 
EXTERNAL  FUNCTN 

C*****  FOR  DOUBLE  PRECISION  CHANGE  REAL  TO  DOUBLE  PRECISION  IN  THE 
C NEXT  STATEMENT . 

DIMENSION  A(NDIM) , B(NDIM) , CENTER(20) 

DIMENSION  WIDTH(20) , WRKSTR( LENWRK) 

INTEGER  DIVAXO,  DIVAXN,  DIVFLG , FUNCLS , IFAIL,  INDEX1 , 

* INDEX2 , J,  K,  LENWRK,  MAXCLS , MAXPTS,  MINPTS,  NDIM, 

* RGNSTR,  RULCLS,  SBRGNS , SBTMPP,  SUBRGN , SUBTMP 
IFAIL=3 

RELERR=1 

FUNCLS=0 

IF(NDIM . LT . 2 . OR . NDIM. GT . 20)  GOTO  300 
IF (MINPTS. GT. MAXPTS)  GOTO  300 


• INITIALISATION  OF  SUBROUTINE 

ZER0=0 

ONE=l 

TWO=2 


HALF=ONE/TWO 
RGNSTR=2*NDIM+3 
ERRMIN  = ZERO 

MAXCLS  = 2**NDIM+2*NDIM**2+6*NDIM+1 

MAXCLS  = MINO (MAXCLS .MAXPTS) 

DIVAXO=0 


C 

c*****  END  SUBROUTINE  INITIALISATION 

IF (MINPTS .LT.O)  SBRGNS=WRKSTR( LENWRK- 1 ) 
IF (MINPTS . LT . 0)  GOTO  280 
DO  30  J=1 , NDIM 

WIDTH(J)=(B(J) -A(J))*HALF 
30  CENTER( J)=A( J )+WIDTH( J ) 

FINEST=ZERO 

WRKS TR( LENWRK )=ZERO 

DIVFLG=1 
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SUBRGN=RGNSTR 

SBRGNS=RGNSTR 

40  CALL  BSRL(NDIM, CENTER, WIDTH, FUNCTN,MAXCLS .RULCLS , 

* ERRMIN , RGNERR , RGNVAL , DIVAXO , DIVAXN) 

FINEST=FINEST+RGNVAL 
WRKSTR ( LENWRK) =WRKSTR ( LENWRK ) +RGNERR 
FUNCLS  = FUNCLS  + RULCLS 
C 

c*****  PLACE  RESULTS  OF  BASIC  RULE  INTO  PARTIALLY  ORDERED  LIST 
C*****  ACCORDING  TO  SUBREGION  ERROR 
IF(DIVFLG.EQ.l)  GO  TO  230 
C 

C*****  WHEN  DIVFLG=0  START  AT  TOP  OF  LIST  AND  MOVE  DOWN  LIST  TREE  TO 
C FIND  CORRECT  POSITION  FOR  RESULTS  FROM  FIRST  HALF  OF  RECENTLY 

C DIVIDED  SUBREGION 

200  SUBTMP=2*SUBRGN 

IF( SUBTMP . GT . SBRGNS ) GO  TO  250 
IF ( SUBTMP . EQ . SBRGNS ) GO  TO  210 
SBTMPP=SUBTMP+RGNSTR 

IF (WRKSTR (SUBTMP) . LT . WRKSTR ( SBTMPP) ) SUBTMP=SBTMPP 
210  IF (RGNERR. GE.WRKSTR( SUBTMP ) ) GO  TO  250 
DO  220  K=1 , RGNSTR 
INDEX1=SUBRGN - K+l 
INDEX2=SUBTMP - K+l 

220  WRKSTR( INDEX1 ) =WRKSTR( INDEX2 ) 

SUBRGN=SUBTMP 
GOTO  200 
C 

C*****  WHEN  DIVFLG=1  START  AT  BOTTOM  RIGHT  BRANCH  AND  MOVE  UP  LIST 
C TREE  TO  FIND  CORRECT  POSITION  FOR  RESULTS  FROM  SECOND  HALF  OF 

C RECENTLY  DIVIDED  SUBREGION 

230  SUBTMP= ( SUBRGN/ (RGNSTR* 2 ) ) *RGNSTR 
IF (SUBTMP.LT. RGNSTR)  GO  TO  250 
IF(RGNERR. LE.WRKSTR( SUBTMP) ) GO  TO  250 
DO  240  K-l, RGNSTR 
INDEX1=SUBRGN-K+1 
INDEX2=SUBTMP - K+l 

240  WRKSTR( INDEX1 )=WRKSTR( INDEX2 ) 

SUBRGN-SUBTMP 
GOTO  230 

C*****  STORE  RESULTS  OF  BASIC  RULE  IN  CORRECT  POSITION  IN  LIST 
250  WRKSTR ( SUBRGN) =RGNERR 
WRKSTR ( SUBRGN - 1 ) =RGNVAL 
WRKSTR ( SUBRGN - 2 ) =DIVAXN 
DO  260  J=1 , NDIM 

SUBTMP=SUBRGN - 2* ( J+l ) 

WRKSTR ( SUBTMP+1 ) =CENTER(J) 

260  WRKSTR ( SUBTMP )=WIDTH(J) 

IF(DIVFLG.EQ.l)  GO  TO  270 

C*****  WHEN  DIVFLG=0  PREPARE  FOR  SECOND  APPLICATION  OF  BASIC  RULE 
CENTER (DIVAXO ) =CENTER (DIVAXO ) +TWO- WIDTH (DIVAXO ) 
SBRGNS=SBRGNS+RGNSTR 
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SUBRGN=SBRGNS 

DIVFLG=1 

CAAAAA  LOOP  BACK  TO  APPLY  BASIC  RULE  TO  OTHER  HALF  OF  SUBREGION 
GO  TO  40 
C 

CAAAAA  END  ORDERING  AND  STORAGE  OF  BASIC  RULE  RESULTS 
C**aaa  MAKE  CHECKS  FOR  POSSIBLE  TERMINATION  OF  ROUTINE 
C 

CAAAAAA  FOR  DOUBLE  PRECISION  CHANGE  ABS  TO  DABS  IN  THE  NEXT  STATEMENT 
270  RELERR=ONE 

IF (WRKSTR (LENWRK)  . LE  . ZERO)  WRKSTR ( LENWRK)  =ZERO 
IF(ABS (FINEST) . NE . ZERO)  RELERR=WRKSTR( LENWRK) /ABS (FINEST) 
IF(RELERR.GT.ONE)  RELERR=ONE 
I F ( S BRGNS+RGNSTR . GT . LENWRK - 2 ) I FAI L=2 
IF ( FUNCLS+FUNCLS “RGNSTR/SBRGNS . GT . MAXPTS ) IFAIL=1 
IF(RELERR . LT . EPS . AND . FUNCLS . GE . MINPTS ) IFAIL=0 
IF(IFAIL.LT. 3)  GOTO  300 
C 

CAAAAA  PREPARE  TO  USE  BASIC  RULE  ON  EACH  HALF  OF  SUBREGION  WITH  LARGEST 
C ERROR 

280  DIVFLG=0 

SUBRGN=RGNSTR 

SUBTMP  = 2*SBRGNS/RGNSTR 

MAXCLS  = MAXPTS /SUBTMP 

ERRMIN  = ABS (FINEST) *EPS /FLOAT ( SUBTMP) 

WRKSTR ( LENWRK) =WRKSTR ( LENWRK) - WRKSTR ( SUBRGN ) 

FINEST=FINEST - WRKSTR ( SUBRGN - 1 ) 

D IVAXO=WRKSTR ( SUBRGN - 2 ) 

DO  290  J-l.NDIM 

SUBTMP=SUBRGN - 2* ( J+l ) 

CENTER ( J ) =WRKSTR ( SUBTMP+1 ) 

290  WIDTH (J ) =WRKSTR ( SUBTMP ) 

WIDTH ( DIVAXO ) =WI DTH ( DIVAXO ) -HALF 

CENTER ( DI VAXO ) =CENTER ( DI VAXO ) -WIDTH (DIVAXO) 

C 

CAAAAA  LOOP  BACK  TO  APPLY  BASIC  RULE 
C 

GOTO  40 
C 

CAAAAA  TERMINATION  POINT 
C 

300  MINPTS=FUNCLS 

WRKSTR ( LENWRK - 1 ) =S  BRGN S 

RETURN 

END 
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SUBROUTINE  BSRL(S , CENTER , HWIDTH , F , MAXVLS , FUNCLS , 

* ERRMIN , ERREST , BASEST , DIVAXO , DIVAXN) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

EXTERNAL  F 
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INTEGER  S,  DIVAXN , DIVAXO,  FUNCLS , INTCLS , I,  MINDEG , MAXDEG , 

* MAXORD,  MINORD , MAXCLS 
DOUBLE  PRECISION  INTVLS 
DIMENSION  CENTER(S) , HUIDTH(S) 

DIMENSION  INTVLS (20) , Z(20),  FULSMS(200),  WEGHTS(200) 

MAXDEG  =12 
MINDEG  = 4 
MINORD  = 0 
ZERO  = 0 
ONE  = 1 
TWO  = 2 
THREE  = 3 
FIVE  = 5 
TEN  = 10 

DO  10  MAXORD  = MINDEG , MAXDEG 

CALL  SYMRL(S , CENTER,  HWIDTH , F,  MINORD,  MAXORD,  INTVLS, 

* INTCLS,  200,  WEGHTS,  FULSMS , IFAIL) 

IF  (IFAIL. EQ. 2)  GOTO  20 

ERREST  = ABS (INTVLS (MAXORD) -INTVLS (MAXORD- 1) ) 

ERRORM  = ABS (INTVLS (MAXORD- 1) -INTVLS (MAXORD- 2) ) 

IF  ( ERREST. NE. ZERO) 

* ERREST  = ERREST*AMAX1 ( ONE/TEN, ERREST/AMAX1 ( ERRE ST/TWO , ERRORM) ) 
IF  (ERRORM. LE.FIVE*ERREST)  GOTO  20 

IF  (2*INTCLS .GT.MAXVLS)  GOTO  20 
IF  ( ERREST. LT.ERRMIN)  GOTO  20 
10  CONTINUE 
20  DIFMAX  = -1 

XI  = ONE/TWO** 2 
X2  = THREE*X1 
DO  30  I = 1,S 
Z(I)  = CENTER ( I ) 

30  CONTINUE 
SUMO  = F(S , Z) 

DO  40  I = 1,S 

Z(I)  = CENTER(I)  - X1*HWIDTH(I) 

SUM1  = F(S,Z) 

Z(I)  = CENTER(I)  + X1*HWIDTH(I) 

SUM1  = SUM1  + F(S,Z) 

Z(I)  = CENTER(I)  - X2*HWIDTH(I) 

SUM2  = F( S , Z) 

Z(I)  = CENTER  ( I ) + X2 -'HWIDTH  ( I ) 

SUM2  = SUM2  + F(S,Z) 

Z(I)  = CENTER(I) 

DIF  = ABS ( ( SUM1 - TWO- SUMO ) - (X1/X2)**2*(SUM2-TWO*SUMO) ) 

IF  (DIF. LT. DIFMAX)  GOTO  40 
DIFMAX  = DIF 
DIVAXN  = I 
40  CONTINUE 

IF  (SUMO . EQ . SUMO+DIFMAX/TWO)  DIVAXN  = MOD(DIVAXO , S)  + 1 
BASEST  = INTVLS (MINORD) 

FUNCLS  = INTCLS  + 4*S 
RETURN 
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END 
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SUBROUTINE  SYMRL(S,  CENTER,  HWIDTH , F,  MINORD , MAXORD , INTVLS , 

* INTCLS , NUMSMS,  WEGHTS,  FULSMS , FAIL) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

MULTIDIMENSIONAL  FULLY  SYMMETRIC  RULE  INTEGRATION  SUBROUTINE 


C 
C 
C 
C 
C 
C 
C 

c kkkkkkkkkkkkkk  PARAMETERS  FOR  SYMRL 
C*****INPUT  PARAMETERS 


THIS  SUBROUTINE  COMPUTES  A SEQUENCE  OF  FULLY  SYMMETRIC  RULE 
APPROXIMATIONS  TO  A FULLY  SYMMETRIC  MULTIPLE  INTEGRAL. 

WRITTEN  BY  A.  GENZ,  MATHEMATICAL  INSTITUTE,  UNIVERSITY  OF  KENT, 
CANTERBURY,  KENT  CT2  7NF , ENGLAND 


k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 


MINORD 


MAXORD 


NUMSMS 


INTEGER  NUMBER  OF  VARIABLES,  MUST  EXCEED  0 BUT  NOT  EXCEED  20 
EXTERNALLY  DECLARED  USER  DEFINED  REAL  FUNCTION  INTEGRAND. 

IT  MUST  HAVE  PARAMETERS  (S,X),  WHERE  X IS  A REAL  ARRAY 
WITH  DIMENSION  S. 

INTEGER  MINIMUM  ORDER  PARAMETER.  ON  ENTRY  MINORD  SPECIFIES 
THE  CURRENT  HIGHEST  ORDER  APPROXIMATION  TO  THE  INTEGRAL, 
AVAILABLE  IN  THE  ARRAY  INTVLS.  FOR  THE  FIRST  CALL  OF  SYMRL 
MINORD  SHOULD  BE  SET  TO  0.  OTHERWISE  A PREVIOUS  CALL  IS 
ASSUMED  THAT  COMPUTED  INTVLS (1),  ...  , INTVLS (MINORD) . 

ON  EXIT  MINORD  IS  SET  TO  MAXORD. 

INTEGER  MAXIMUM  ORDER  PARAMETER,  MUST  BE  GREATER  THAN  MINORD 
AND  NOT  EXCEED  20.  THE  SUBROUTINE  COMPUTES  INTVLS (MINORD+1) , 
INTVLS (MINORD+2) , . . . , INTVLS (MAXORD) . 

REAL  ARRAY  OF  DIMENS ION (MAXORD)  OF  GENERATORS. 

ALL  GENERATORS  MUST  BE  DISTINCT  AND  NONNEGATIVE. 

INTEGER  LENGTH  OF  ARRAY  FULSMS,  MUST  BE  AT  LEAST  THE  SUM  OF 
THE  NUMBER  OF  DISTINCT  PARTITIONS  OF  LENGTH  AT  MOST  S 
OF  THE  INTEGERS  0,1,..., MAXORD - 1 . AN  UPPER  BOUND  FOR  NUMSMS 
WHEN  S+MAXORD  IS  LESS  THAN  19  IS  200 
C******OUTPUT  PARAMETERS 

C INTVLS  REAL  ARRAY  OF  D I MEN SI ON (MAXORD) . UPON  SUCCESSFUL  EXIT 

INTVLS (1) , INTVLS (2) , . . . , INTVLS (MAXORD)  ARE  APPROXIMATIONS 
TO  THE  INTEGRAL.  INTVLS (D+l)  WILL  BE  AN  APPROXIMATION  OF 
POLYNOMIAL  DEGREE  2D+1. 

INTEGER  TOTAL  NUMBER  OF  F VALUES  NEEDED  FOR  INTVLS (MAXORD) 
REAL  WORKING  STORAGE  ARRAY  WITH  DIMENSION  (NUMSMS).  ON  EXIT 
WEGHTS (J)  CONTAINS  THE  WEIGHT  FOR  FULSMS (J). 

REAL  WORKING  STORAGE  ARRAY  WITH  DIMENSION  (NUMSMS).  ON  EXIT 
FULSMS (J)  CONTAINS  THE  FULLY  SYMMETRIC  BASIC  RULE  SUM 
INDEXED  BY  THE  JTH  S- PARTITION  OF  THE  INTEGERS 
0,1, . . . , MAXORD- 1. 

INTEGER  FAILURE  OUTPUT  PARAMETER 

FAIL=0  FOR  SUCCESSFUL  TERMINATION  OF  THE  SUBROUTINE 
FAIL=1  WHEN  NUMSMS  IS  TOO  SMALL  FOR  THE  SUBROUTINE  TO 

CONTINUE.  IN  THIS  CASE  WEGHTS(l),  WEGHTS(2),  ..., 
WEGHTS ( NUMSMS ) , FULSMS ( 1 ) , FULSMS ( 2 ) , 


INTCLS 

WEGHTS 

FULSMS 


FAIL 


40 


c 

c 

c 

c 

c 


FULSMS ( NUMSMS ) AND  INTVLS(l),  INTVLS ( 2 ) , . . . , 
INTVLS(J)  ARE  RETURNED,  WHERE  J IS  MAXIMUM  VALUE  OF 
MAXORD  COMPATIBLE  WITH  THE  GIVEN  VALUE  OF  NUMSMS. 


FAIL=2  WHEN  PARAMETERS  S.MINORD,  MAXORD  OR  G ARE  OUT  OF 
RANGE 


EXTERNAL  F 

C***  FOR  DOUBLE  PRECISION  CHANGE  REAL  TO  DOUBLE  PRECISION 


C 


IN  THE  NEXT  STATEMENT 


INTEGER  D,  I,  FAIL,  K(20) , INTCLS , PRTCNT , L,  M(20),  MAXORD, 
* MINORD , MODOFM,  NUMSMS,  S,  SUMCLS 
DOUBLE  PRECISION  INTVLS,  INTMPA,  INTMPB 
DOUBLE  PRECISION  MOMTOL,  INTVAL , MOMENT,  MOMPRD , MOMNKN 
DIMENSION  INTVLS (MAXORD) , CENTER(S),  HWIDTH(S) 

DIMENSION  FULSMS (NUMSMS ) , WEGHTS (NUMSMS ) 

DIMENSION  MOMPRD(20 , 20) , MOMENT ( 20) , G(20) 

C PATTERSON  GENERATORS 

DATA  G ( 1 ) , G(2)  /0 . 0000000000000000 , 0..  7745966692414833/ 

DATA  G(3) , G(4)  /0 . 9604912687080202 , 0 . 4342437493468025/ 

DATA  G(5) , G ( 6 ) /0 . 9938319632127549 , 0 . 8884592328722569/ 

DATA  G(7) , G ( 8 ) /0 . 6211029467372263 , 0 . 2233866864289668/ 

DATA  G(9) , G(10) , G(ll) , G(12)  /0 . 1 , 0.2,  0.3,  0.4/ 

C 

C***  PARAMETER  CHECKING  AND  INITIALISATION 
FAIL  = 2 
MAXRDM  =20 
MAXS  = 20 

IF  (S . GT .MAXS  .OR.  S.LT.l)  RETURN 
IF  (MINORD. LT.O  .OR.  MINORD . GE .MAXORD)  RETURN 
IF  ( MAXORD. GT. MAXRDM)  RETURN 
ZERO  = 0 
ONE  = 1 
TWO  = 2 
MOMTOL  = ONE 
10  MOMTOL  = MOMTOL/TWO 

IF  (MOMTOL+ONE . GT . ONE)  GO  TO  10 
HUNDRD  = 100 

MOMTOL  = HUNDRD*TWO*MOMTOL 
D = MINORD 

IF  (D.EQ.O)  INTCLS  = 0 

C***  CALCULATE  MOMENTS  AND  MODIFIED  MOMENTS 
DO  20  L=l, MAXORD 
FLOATL  = L + L - 1 
MOMENT (L)  = TWO/FLOATL 
20  CONTINUE 

IF  (MAXORD. EQ.l)  GO  TO  50 
DO  40  L=2, MAXORD 

INTMPA  = MOMENT (L-l) 

GLSQRD  = G(L-1)**2 
DO  30  I=L, MAXORD 
INTMPB  = MOMENT(I) 

MOMENT ( I ) = MOMENT ( I ) - GLSQRD*INTMPA 
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INTMPA  = INTMPB 
30  CONTINUE 

IF  (MOMENT (L)**2.LT. (MOMTOL*MOMENT(l) )**2)  MOMENT (L)  = ZERO 
40  CONTINUE 
50  DO  70  L=1 , MAXORD 

IF  (G(L) . LT . ZERO)  RETURN 

MOMNKN  = ONE 

MOMPRD(L, 1)  = MOMENT ( 1 ) 

IF  (MAXORD. EQ„ 1)  GO  TO  70 
GLSQRD  = G(L)**2 
DO  60  1=2 .MAXORD 

IF  (I.LE.L)  GISQRD  = G(I-1)**2 
IF  (I.GT.L)  GISQRD  = G(I)**2 
IF  (GLSQRD. EQ. GISQRD)  RETURN 
MOMNKN  = MOMNKN/ (GLSQRD -GISQRD) 

MOMPRD(L.I)  = MOMNKN-MOMENT ( I ) 

60  CONTINUE 
70  CONTINUE 
FAIL  = 1 
C 

C***  BEGIN  LOOP  FOR  EACH  D 

C FOR  EACH  D FIND  ALL  DISTINCT  PARTITIONS  M WITH  MOD(M))=D 

C 

80  PRTCNT  = 0 
INTVAL  = ZERO 
MODOFM  = 0 

CALL  NXPRT( PRTCNT,  S,  M) 

90  IF  ( PRTCNT. GT.NUMSMS)  RETURN 
C 

C***  CALCULATE  WEIGHT  FOR  PARTITION  M AND  FULLY  SYMMETRIC  SUMS 
C***  WHEN  NECESSARY 
C 

IF  (D.EQ. MODOFM)  WEGHTS ( PRTCNT)  = ZERO 
IF  (D.EQ. MODOFM)  FULSMS ( PRTCNT)  = ZERO 
FULWGT  = WHT ( S , MOMENT , M , K , MODOFM , D , MAXRDM , MOMPRD ) 

SUMCLS  = 0 

IF  (WEGHTS (PRTCNT) .EQ. ZERO  .AND.  FULWGT . NE . ZERO)  FULSMS (PRTCNT) 
* FLSM(S , CENTER,  HWIDTH , MOMENT,  M,  K,  MAXORD,  G,  F,  SUMCLS) 
INTCLS  = INTCLS  + SUMCLS 
INTVAL  = INTVAL  + FULWGT*FULSMS( PRTCNT) 

WEGHTS (PRTCNT)  = WEGHTS ( PRTCNT)  + FULWGT 
CALL  NXPRT( PRTCNT,  S,  M) 

IF  (M(l) .GT. MODOFM)  MODOFM  = MODOFM  + 1 
IF  (MODOFM. LE.D)  GO  TO  90 
C 

C***  END  LOOP  FOR  EACH  D 

IF  (D.GT.O)  INTVAL  = INTVLS(D)  + INTVAL 
INTVLS (D+l)  = INTVAL 
D = D + 1 

IF  (D.LT. MAXORD)  GO  TO  80 
C 

C***  SET  FAILURE  PARAMETER  AND  RETURN 
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FAIL  = 0 
MINORD  = MAXORD 
RETURN 
END 

SUBROUTINE  NXPRT(PRTCNT , S,  M) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C 

C***  SUBROUTINE  TO  COMPUTE  THE  NEXT  S PARTITION 
C 


INTEGER  S,  M(S) , PRTCNT , I,  MSUM 
IF  (PRTCNT. GT.O)  GO  TO  20 
DO  10  1=1, S 
M( I ) = 0 
10  CONTINUE 
PRTCNT  = 1 
RETURN 

20  PRTCNT  = PRTCNT  + 1 
MSUM  = M(l) 

IF  (S.EQ.l)  GO  TO  60 
DO  50  1=2 ,S 

MSUM  = MSUM  + M(I) 

IF  (M( 1) . LE . M(I )+l)  GO  TO  40 
M(l)  = MSUM  - (I-1)*(M(I)+1) 

DO  30  L=2 , I 

M(L)  = M(I ) + 1 
30  CONTINUE 
RETURN 

40  M( I ) = 0 

50  CONTINUE 
60  M(l)  = MSUM  + 1 
RETURN 
END 
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DOUBLE  PRECISION  FUNCTION  WHT 

(S,  INTRPS , M,  K,  MODOFM,  D,  MAXRDM , MOMPRD) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C***  SUBROUTINE  TO  CALCULATE  WEIGHT  FOR  PARTITION  M 
C 

INTEGER  S,  M(S) , K(S),  D,  MAXRDM,  MI,  KI , Ml,  K1 , MODOFM 

DOUBLE  PRECISION  INTRPS,  MOMPRD 

DIMENSION  INTRPS ( S ) , MOMPRD ( MAXRDM , MAXRDM ) 

ZERO  = 0 
DO  10  1=1, S 

INTRPS (I)  = ZERO 
K(I)  = 0 
10  CONTINUE 

Ml  = M( 1 ) + 1 
Kl  = D - MODOFM  + Ml 
20  INTRPS (1)  = MOMPRD (Ml, Kl) 

IF  (S.EQ.l)  GO  TO  40 


kkk 
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DO  30  1=2, S 
MI  = M( I ) + 1 
KI  = K( I ) + MI 

INTRPS ( I ) = INTRPS(I)  + MOMPRD(MI , KI)*INTRPS ( I - 1) 

INTRPS(I-l)  = ZERO 
Kl  = Kl  - 1 
K(I)  = K( I ) + 1 
IF  (K1.GE.M1)  GO  TO  20 
Kl  = Kl  + K( I) 

K( I ) = 0 
30  CONTINUE 
40  WHT  = INTRPS (S) 

RETURN 

END 
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DOUBLE  PRECISION  FUNCTION  FLSM 

( S , CENTER , HWIDTH , X , M , MP , MAXORD , G , F, SUMCLS ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

EXTERNAL  F 
C 

CAAA  FUNCTION  TO  COMPUTE  FULLY  SYMMETRIC  BASIC  RULE  SUM 
C 

INTEGER  S,  M(S) , MP(S) , MAXORD,  SUMCLS,  IXCHNG , LXCHNG , I,  L, 

* IHALF , MPI , MPL 
DOUBLE  PRECISION  INTWGT,  INTSUM 
DIMENSION  G( MAXORD) , X(S) 

DIMENSION  CENTER(S) , HWIDTH(S) 

ZERO  = 0 
ONE  = 1 
TWO  = 2 
INTWGT  = ONE 
DO  10  1=1, S 
MP ( I ) = M( I ) 

IF  (M(I).NE.O)  INTWGT  = INTWGT/TWO 
INTWGT  = INTWGT*HWIDTH(I) 

10  CONTINUE 
SUMCLS  = 0 
FLSM  = ZERO 
C 

C A A A A A A A COMPUTE  CENTRALLY  SYMMETRIC  SUM  FOR  PERMUTATION  MP 
20  INTSUM  = ZERO 
DO  30  1=1, S 

MPI  = MP ( I ) + 1 

X(I)  = CENTER(I)  + G(MPI)*HWIDTH(I) 

30  CONTINUE 

40  SUMCLS  = SUMCLS  + 1 

INTSUM  = INTSUM  + F(S,X) 

DO  50  1=1, S 

MPI  = MP ( I ) + 1 

IF(G(MPI) .NE.ZERO)  HWIDTH(I)  = -HWIDTH(I) 
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X(I)  = CENTER(I)  + G(MPI)*HWIDTH(I) 

IF  (X ( I ) . LT . CENTER ( I ) ) GO  TO  40 
50  CONTINUE 

C*******  END  INTEGRATION  LOOP  FOR  MP 
C 

FLSM  = FLSM  + INTWGT*INTSUM 
IF  (S.EQ.l)  RETURN 
C 

C*******  FIND  NEXT  DISTINCT  PERMUTATION  OF  M AND  LOOP  BACK 
C TO  COMPUTE  NEXT  CENTRALLY  SYMMETRIC  SUM 

DO  80  1=2, S 

IF  (MP(I-l) .LE.MP(I))  GO  TO  80 
MPI  = MP ( I ) 

IXCHNG  =1-1 
IF  (I.EQ.2)  GO  TO  70 
I HALF  = IXCHNG/2 
DO  60  L=1 , IHALF 
MPL  = MP(L) 

IMNUSL  = I - L 
MP(L)  = MP ( IMNUSL) 

MP ( IMNUSL)  = MPL 

IF  (MPL. LE. MPI)  IXCHNG  = IXCHNG  - 1 
IF  (MP (L) . GT . MPI ) LXCHNG  = L 
60  CONTINUE 

IF  (MP ( IXCHNG ) .LE. MPI)  IXCHNG  = LXCHNG 
70  MP ( I ) = MP( IXCHNG) 

MP( IXCHNG)  = MPI 
GO  TO  20 
80  CONTINUE 

c*****  END  LOOP  FOR  PERMUTATIONS  OF  M AND  ASSOCIATED  SUMS 
C 

RETURN 

END 

C ~k  ~k  ~k  ~k  ~k  -k  -k  'k  V?  -k  ~k  -k  ~k  -k  k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 


C k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k -.V  k k -V  k k k k k k k k k k k k k k k k k 'A' 


c*** 

c*** 


DOUBLE  PRECISION  FUNCTION  QFCLG(NDIM, Z) 

This  function  computes  the  convective  heat  transfer  flux  to 
the  ceiling  at  location  (X, Y)=(Z(1) , Z(2) ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DIMENSION  Z(2) 

COMMON/AINTCH/H , HTCT , THT , THTQHP , Cl , C2 , C3 , XF , YF , TC 
X-Z(l) 

Y=Z(2) 

R= ( (X-XF) **2 . D0+ ( Y- YF) **2 . DO ) ** ( 1 . DO/2 . DO ) 

RDH=R/H 

FF=( 1 . DO  - 1 . 1D0*(RDH**0 . 8 DO ) +0 . 808D0* (RDH**1 . 6 DO ) ) / 

+ (1 . DO-1 . 1D0*(RDH**0 . 8D0)+2 . 2D0*(RDH**1 . 6D0)+ 

+ 0.69D0*(RDH**2.4D0)) 

IF(RDH . LT . 0 . 2D0)THEN 

HTCLDH=C1* ( 1 . DO  - C2*RDH) 

TADDIM=10 . 22D0 - 14 . 9D0*RDH 
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ELSE 


HTCLDH-C3* ( 1 . DO/ (RDH**1 . 2D0 ) ) * (RDH- 0 . 077 IDO ) / 

+ (RDH+O . 279D0) 

TADDIM=8 . 390913361D0*FF 

ENDIF 

HTCL=HTCLDH*HTCT 
TAD=TADDIM*THTQHP+THT 
QFCLG-HTCL* ( TAD - TC ) 

RETURN 

END 

Q 'k  •Jt  ~k  ~,V  ~k  ~k  ~k  'k  k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k 
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Figure  1.  Fire  in  a rectangular  parallelopiped  enclosure. 
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Figure  2.  The  fire  and  the  equivalent  source  in  the  lower  layer  and  the  continuation  source  in  the 
extended  upper  layer. 
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Figure  3.  Overview  of  the  CEILIIT  Algorithm. 
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Figure  4.  A plot  of  dimensionless  ceiling  jet  velocity  distribution,  VCJ/V MAX’  as  a functi°n  °f 
zj (0.235)  per  Eq.  (29). 
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Figure  5.  Plots  of  dimensionless  ceiling  jet  temperature  distribution,  0,  as  a function  of  z/(0.23S) 
per  Eq.  (32)  for  cases  when  ©s  is  < 0,  between  0 and  1,  and  > 0. 
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Figure  6. 


Definition  of  the  Wall  Stagnation  Points,  Room  Corners,  and  Wall  Segments. 
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